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1.0  Introduction 

This  report  documents  a  program  for  the  solution  of  algebraic  and 
implicitly  defined  stiff  differential  equations.   We  were  particularly 
interested  in  solving  very  large  systems  of  differential  equations 
arising  from  partial  differential  equations  where  the  finite  element 
method  has  been  applied  in  the  space  variables. 

Our  original  goal  was  to  use  a  compact  storage  scheme  for  the 
large  matrices  involved  and  to  use  iteration  to  solve  the  linear 
algebraic  systems  which  occur.   However,  the  resulting  program  is  easily 
adapted  to  different  applications  through  user  modifications  accomplished 
by  replacement  of  one  or  two  relatively  simple  subroutines.   Thus  the 
program  is  a  powerful  one  which  can  be  used  in  a  variety  of  applications. 
Four  examples,  illustrating  different  matrix  storage  techniques  and 
different  linear  equation  solvers  are  given  in  the  appendix.   Other 
storage  schemes  and  solution  methods,  e.g.  symmetric  Jacobian  with 
Cholesky  decomposition,  are  relatively  simple  to  implement. 

In  Section  2  a  brief  review  of  the  integration  scheme  is  given. 
In  Section  3  a  discussion  of  differences  between  this  program  and  the 
one  from  which  it  was  adapted  is  given.   Section  4  is  devoted  to  a 
discussion  of  information  concerning  the  use  of  this  package. 

2.0  Theoretical  Background 

Consider  the  system  of  N  =  m  +  I     differential  and  algebraic 
equations  in  y± ym»vi»  *  *  *  »  V£» 

(1)  F(y,y,t)  +  P(t)  V  =  0  , 


with  all  or  some  of  the  initial  values  y  (t0),...,y  (t~) ,v  (t„) , . . .  , 
v  (t~)   specified.   Enough  of  the  above  values  must  be  given  in  order 
to  determine  the  remaining  values  and  initial  values  for  any  of  the 
derivatives,   y-i>...,y   which  appears  in  equation  (1).   In  equation  (1) , 
P(t)   is  an  N  x  I     matrix,   F  is  a  vector  with  N  components,  and 

V  is  a  vector. 

The  program  documented  here  is  a  modification  of  a  program  due  to 
Brown  and  Gear  [1],   The  method  of  solution  is  a  modification  of  Gear's 
method  for  stiff  differential  equations  [2],   The  application  to  differ- 
ential algebraic  systems  was  given  by  Gear  [3],   We  will  briefly  describe 
the  method  here  for  completeness,  and  refer  the  reader  to  the  references 
for  more  details. 

Suppose  that  approximate  solution  values  are  known  at  a  number  of 

equally  spaced  points,   t  _ ,t  _~ t  ,  ,   and  are  represented  by 

y     ,...,y       respectively.   Let  V(t  -)   be  represented  by 

V  .   Use  of  a  backward  differentiation  formula  gives 

,  «(n)    1   ,    (n)       (n-1)  (n-k)  . 

hy    "J  (aQ  y    +  a-L  y    '•+■...+  afc  y      )  , 

where  h  =  t.in  -  t.  .   The  coefficients  a.   and  3--,  are  from  Gear 
l+l    l  i       0 

[2]  ,  p.  217.   Substitution  of  this  into  (1)  gives 

(n)  a0  (n)    ,    r  i   n   _l  T>/fc   n   ttC11)    _  n 


(2)  F(yW,    --SL.     yW    +  I      ,    t   )   +p(t   ) 

J  8nh     J  Ln   *     n  n 


Bo 


V 


as   the  equation  which     y  and     V  must   satisfy.      In  equation    (2), 


v  1      t  (n-1)  (n-k). 

^n  =  6Th   (al  y  +  •  •  •  +  ak  y  )    . 


In  general,  equation  (2)  represents  a  system  of  nonlinear  equations  for 

;ebrai( 
(n),0 


y     and  V    .   The  method  used  for  solving  this  system  of  algebraic 


equations  is  a  variant  of  Newton's  method.   The  initial  guess,   y 

is  obtained  by  polynomial  extrapolation  using  a  Hermite  polynomial 
interpolating  the  known  values  y     ,  y     ,  ...,  y      .   Thus 
the  predicted  values  has  the  form 


._.    (n),0   .  -   -(n-1)   -   (n-1)         -     (n-k) 

(3)   y   '  =hS,  y      +  a,  y'     +  ...  +  a  .  yK 
±  1  n-K 


The  application  of  Newton's  Method  to  equation  (2)  then  yields  the 
corrector  equation 

(n),i+l  _   (n),i 

)=-F^(n),i'  -Fh^'1-  V^V  v(n),i  • 

r(n),i+l   _<n)ti  /  V 


V 


-  V 


where  J  is  the  Jacobian  matrix, 


(5) 


3y   6Qh  9y 


Gear  shows  that  the  initial  guess  for  V     is  not  important , 
and  thus  V       is  used.   Up  to  three  iterations  are  performed  on  the 
corrector  equation.   The  matrix  J  is  not  evaluated  at  each  iteration, 
nor  even  at  each  timestep.   J  is  evaluated  whenever  (i)  the  timestep 
or  order  is  changed,  or  (ii)  the  corrector  iteration  fails  to  converge 
in  three  iterations. 

If  the  corrector  iteration  fails  to  converge,  the  J  matrix  is 
evaluated,  unless  it  had  already  been  evaluated  at  the  current  time. 


If  it  has  been  evaluated  at  the  current  time,  the  timestep  is  reduced 
by  a  factor  of  4.   In  either  case  the  step  is  then  retired. 

If  the  corrector  iteration  converges,  the  local  error  is  estimated, 
based  on  the  fact  that  the  local  error  is  approximately  proportional  to 
the  difference  between  the  predicted  and  corrected  values  of  y    . 
For  this  purpose  a  relative  error  tolerance  is  used  for  large  solutions 
and  an  absolute  error  tolerance  for  small  solutions.   The  root-mean- 
square  norm  (Euclidean  norm  divided  by  the  square  root  of  the  number  of 
components)  is  used  for  the  vector  with  components  e./ymax.  ,  where 

e.   is  the  estimated  local  error  in  y.    and  ymax.  =  max  (|y.   | ,  1)  . 
1  X  1  0<k<n   1 

If  the  error  is  larger  than  that  specified  by  the  user,  an  acceptable 
timestep  is  estimated  for  the  current  order  or  order  one  lower,  and 
the  step  repeated.   Up  to  three  such  failures  are  permitted,  after  which 
an  attempt  is  made  to  start  over  with  a  first  order  method. 

When  using  a  method  of  order  q  ,  the  program  takes  at  least 
q  +  1  steps  before  changing  the  timestep.   Changes  in  timestep  are 
preceeded  by  calculation  of  the  predicted  timestep  at  current  order  and 
order  one  higher  and  one  lower.   If  the  timestep  can  be  increased  by  more 
than  10%,  the  order  corresponding  to  the  largest  possible  timestep  is 
used.   If  the  timestep  cannot  be  increased  by  at  least  10%,  the  current 
order  and  stepsize  are  retained  for  at  least  10  more  steps. 

After  each  step  a  test  is  made  to  determine  whether  time  has 
advanced  to  or  beyond  the  input  end  time.   Control  is  returned  to  the 
calling  program  if  it  has. 


At  the  initial  call,  no  history  of  the  solution  is  available,  so 
the  program  must  begin  with  a  first  order  method,  taking  two  such  steps 
in  accordance  with  the  above  description.   The  timestep  must  be  suitably 
small,  again  in  accordance  with  the  above.   At  the  point  the  program  can 
begin  to  increase  the  order  of  the  method  and  the  timestep.   Because 
the  Jacobian  matrix  must  be  generated  whenever  the  timestep  is  changed, 
it  is  not  efficient  to  try  too  large  a  timestep  initially.   Because  the 
program  rapidly  finds  the  best  order  and  timestep,  it  is  relatively  cheap 
to  underestimate  the  initial  timestep  compared  to  the  cost  of  overestimating 
it. 

3.0  Differences  compared  with  DFASUB 

The  principal  differences  between  the  SDESOL/LDASUB  package  and 
DFASUB,  and  the  reasons  for  incorporating  them  are  as  follows. 
(i)   The  main  goal  of  this  revision  was  to  generate  a  program  which 
could  solve  very  large  sparse  systems  of  differential  equations  efficiently, 
both  in  terms  of  storage  requirements  and  execution  time.   We  are 
particularly  interested  in  the  solution  of  ordinary  differential  equa- 
tions arising  from  partial  differential  equations  where  the  finite 
element  method  has  been  used  to  discretize  the  problem  in  space. 

Large  sparse  problems  require  at  least  a  different  system  of  stor- 
ing the  Jacobian  matrix  and  possibly  the  use  of  iteration  to  solve  for 
the  quasi-Newton  iterates  in  equation  (2.4).   Two  such  subroutine 
packages,  to  be  used  with  the  basic  subroutines,  have  been  provided. 
Another  package  using  standard  elimination  techniques  is  also  provided 
and  is  convenient  for  smaller  systems  of  equations.   Use  of  any  of  these 


options  requires  the  user  to  supply  a  subroutine  to  evaluate  his  form 
of  the  equations  (1),  and  for  efficiency,  a  subroutine  to  explicitly 
evaluate  the  Jacobian.   A  subroutine  to  approximate  a  full  Jacobian 
through  numerical  differencing  has  been  provided.   With  the  exception 
of  a  minor  correction,  this  is  the  same  as  given  in  [1],   While  use  of 
this  routine  is  convenient,  it  is  inefficient  and  should  be  avoided 
for  large  systems.   It  is  anticipated  that  the  user  can  provide  his  own 
subroutine  package,  using  his  own  storage  scheme  for  the  Jacobian,  and 
with  a  suitable  equation  solver  for  the  Newton  iterates.   There  should 
be  no  need  to  disturb  the  basic  package  which  carries  out  the  time 
integration. 

(ii)   For  user  convenience,  without  a  major  rewrite  of  DFASUB,  a  driver 
routine,  SDESOL,  to  be  referenced  by  the  user  and  which  then  communicates 
with  LDASUB  was  written.   The  chief  function  of  SDESOL  is  to  set  up  a 
number  of  references  to  work  storage  areas  required  by  LDASUB.   In 
addition,  some  testing  of  parameters  is  accomplished,  and  a  subroutine 
to  calculate  initial  values  of  derivatives  is  called, 
(iii)   A  subroutine  to  calculate  initial  values  of  derivatives,  DERVAL, 

has  been  provided.   The  routine  provided  requires  that  the  first  m 

8F 
rows  of  —^    be  nonsingular,  which  does  not  need  to  be  true  in  the 

general  case.   For  this  reason,  and  others  discussed  in  Section  4,  the 

user  may  need  to  provide  either  his  own  version  of  DERVAL  to  evaluate 

the  derivatives  initially,  or  else  he  may  supply  initial  values  and  a 

dummy  version  of  DERVAL. 

(iv)   Other  changes  made  in  generating  LDASUB  from  DFASUB  were  to  simplify 

the  code  for  the  particular  type  of  problem  we  wish  to  solve,  while 


others  were  to  enhance  the  usability  of  the  code.   Some  errors  were 
also  corrected,  notably  two  errors  in  coefficients  for  the  fifth  and 
sixth  order  methods.   DFASUB  had  the  capability  of  computing  various 
elements  of  the  Jacobian  at  different  times  if  they  had  different 
dependencies,  with  the  possibility  of  inverting  that  part  of  the  matrix 
at  that  time,  if  it  could  be  done.   This  could  result  in  increased 
efficiency  in  certain  problems,  at  some  expense  in  convenience,  but  for 
our  purposes  it  was  not  considered  useful,  and  was  removed.   Therefore 
only  one  call  is  made  to  evaluate  the  Jacobian.   The  Jacobian  was 
evaluated  at  the  beginning  of  each  timestep  in  DFASUB,  but  this  has 
been  eliminated  in  LDASUB,  in  accordance  with  the  description  in  Section 
2.   A  subroutine,  S2,  was  called  from  DFASUB  to  evaluate  time  dependent 
terms  whenever  time  was  changed.   This  is  reasonable,  since  the  function 
evaluation  routine  may  be  called  several  times  at  a  given  value  of  the 
independent  variable.  We  have  removed  this,  preferring  to  test  for  a 
new  time  in  the  routines  where  time  dependent  terms  appear,  then 
evaluating  and  storing  them  internally  to  that  routine  when  necessary. 
This  helps  make  the  function  evaluation  more  self  contained,  as  well. 

In  DFASUB  extra  parameters  in  the  calling  sequence  allowed  the 
user  to  communicate  constants  to  the  function  evaluation  and  Jacobian 
subroutines  via  DFASUB.   We  believe  this  is  inefficient  and  confusing, 
and  removed  this  capability,  preferring  to  communicate  from  the  main 
program  to  these  subroutines  via  Common,  or  possibly  through  multiple 
entry  points. 

The  norm  used  for  error  tests  in  DFASUB  is  the  Euclidean  norm. 
This  has  the  undesirable  property  that  for  large  systems  the  allowable 


error  criterion  may  be  large.   We  therefore  changed  to  the  root-mean- 
square  norm  in  LDASUB,  which  is  simply  the  Euclidean  norm  divided  by 
the  square  root  of  the  number  of  components.   One  other  change  was 
made  in  the  error  tests.   As  noted  in  Section  2,  the  error  vector  has 
components  e./ymax.  ,  where  e.   is  the  estimated  local  error  in  the 

i —  variable   y.  ,  and  ymax,  =  max  (|y.   |,  1)  .   In  DFASUB  the 
1  0 <k<n    1 

maximum  was  taken  only  up  to  the  previous  timestep,  n-1  .   This  change 

was  incorporated  because  some  of  the  problems  in  which  we  were  interested 

began  with  many  components  at  zero,  but  which  very  rapidly  became  large, 

12 
around  10    or  more.  Without  updating  the  value  of  ymax.  ,  the 

size  of  the   timestep  was  artificially  kept  extremely  small  in  order  to 
satisfy  an  unreasonable  error  tolerance.   For  this  reason,  the  maximum 
value  of  the  component  was  updated  before  the  norm  of  the  relative  errors 
was  computed.   For  problems  where  values  range  near  to  one,  the  modifica- 
tion will  result  in  no  appreciable  change  in  performance. 

The  printout  of  counters,  timestep,  time,  and  values  of  the  dependent 
variables  was  made  an  option  through  a  parameter  in  the  calling  sequence. 
An  additional  value  printed  is  the  order  of  the  method  being  used. 

DFASUB  incorporated  the  capability  of  terminating  if  a  certain 
number  of  floating  point  overflows  had  occurred.   This  capability  was 
removed  from  LDASUB. 

The  final  modification  to  the  program  was  the  incorporation  of  a 
restart  capability  without  having  to  begin  again  with  a  first  order 
method.   This  was  accomplished  by  adding  two  entry  points  to  LDASUB. 
One,  LDASAV,  saves  values  internal  to  LDASUB,  returning  them  to  the 
main  program,  where  they  can  be  saved  for  the  time  at  which  the  calculation 


is  to  be  resumed.   At  that  time,  another  entry  point,  LDARST,  restores 
those  values  internal  to  LDASUB,  while  other  necessary  values  are 
restored  through  the  calling  sequence. 

4.0  Subroutine  Descriptions 

The  description  of  subroutines  is  divided  into  two  subsections. 
The  first  deals  with  the  basic  integration  routine  and  other  subroutines 
which  make  up  the  core  package,  and  which  should  not  be  modified  by  the 
casual  user.   The  second  deals  with  a  set  of  supporting  subroutines,  at 
least  one  of  which  must  be  supplied  by  the  user  since  it  defines  his 
system  of  equations.   The  others  may  be  usable  in  the  form  given  in 
one  of  the  examples,  or  can  be  defined  by  the  user  to  accomplish  his 
desired  implementation. 

4.1  Basic  Subroutine  Package 

4.1.1   Subroutine  SDESOL 

This  routine  is  the  only  one  which  needs  to  be  referenced  by  most 
users.   It  serves  as  a  driver  for  the  integration  routine,  LDASUB. 
SDESOL  has  a  simpler  calling  sequence  than  LDASUB  and  relieves  the 
user  of  having  to  set  up  a  number  of  auxiliary  storage  arrays.   In 
addition,  the  routine  calls  DERVAL  to  calculate  the  values  of  the 
derivatives  on  the  initial  call. 

The  calling  sequence  is 
CALL   SDESOL (Y ,YL , T , TEND, NY ,NL ,M, JSKF ,MAXDER, IPRT ,H, HMIN ,HMAX , RMSEPS , W) 
where  the  parameters  are  defined  as  follows. 


Y  -  Input  and  output.   An  array  dimensioned  (7, NY).   On  the  initial 
call  this  array  contains  the  initial  values  of  the  dependent 

variables  v.,  i=l ,  . . . ,  m  in  Y(l,i)  .   During  execution 

d3yi   hj 
of  the  program  the  approximate  values  of  r—  •  -rr     is 

d  t J    J  * 
stored  in  Y(j+l,i)  .   Here  h  is  the  current  stepsize. 

These  values  must  not  be  changed  between  returns  to  the 

calling  program  and  subsequent  entries  to  SDESOL.   It  is 

possible  to  interpolate  for  values  of  the  dependent  variables 

at  times  other  than  those  calculated  by  using  the  formula 

y.(t+s)  =  I     Y(j+l,i)  (  t-  )   ,  where  q   is  the  order  formula 
1       j=0        V n  / 

currently  in  use,  and  is  obtained  as  q  =  |JSKF/10|  . 

YL  -  input  and  output.   Array  of  linear  variables,  v.,  i=l,  ...,  I   . 
The  user  supplies  initial  values  for  these  variables,  and 
during  execution  it  contains  current  values  of  the  linear 
variables. 

T  -  input  and  output.   The  user  supplies  the  initial  time,  which 
is  updated  to  current  time  during  execution. 

TEND  -input.  Time  at  which  the  integration  is  to  end.   This  is  the 
only  parameter  normally  changed  by  the  user  between  succes- 
sive entries  to  SDESOL. 

NY  -  input.   The  number  of  differential  and  nonlinear  variables,  m  . 

NL  -  input.   The  number  of  linear  variables,   £  .   This  may  be  zero. 

M  -   input.   The  number  of  variables  to  be  included  in  the  local 

error  test.   The  error  test  will  be  performed  for  the  variables 
y. ,  i=l,  ...5  M  .   The  M  used  is  no  greater  than  NY  . 
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JSKF  -input  and  output.   An  indicator:  on  input, 

JSKF  =  0   indicates  that  this  is  the  initial  call  to  SDESOL. 

Initial  values  of  the  derivatives  are  calculated  and  auxiliary 

storage  references  are  set  up.   This  also  indicates  to 

subroutine  LDASUB  to  initialize  parameters  and  begin  with 

a  first  order  integration  method. 

JSKF  >  0   indicates  a  continuation  of  a  previous  call  to  SDESOL 

JSKF  =  -  1   indicates  a  restart  call.   This  is  discussed 

further  in  Section  4.1.2. 

JSKF  <  -  1  may  result  from  the  user  neglecting  to  test  for 

error  returns  from  SDESOL.   Because  of  this  possibility, 

the  run  is  terminated  with  an  appropriate  comment  when 

JSKF  <  -  1  is  input. 

On  output,  JSKF  normally  is  a  two  digit  number,   ±  qp  .   q 

is  the  order  of  the  formula  currently  being  used  for  the 

integration,   p  is  an  indicator  determining  the  type  of 

return.   JSKF  >  0  ,  p  =  1  is  the  normal  return.   Note  that 

SDESOL  may  be  re-entered  to  continue  the  solution  without 

changing  JSKF.   JSKF  <  0  is  an  error  return,  with  the  value 

of  p  indicating  the  error,  as  follows. 

p  =  1  error  test  failure  for  H  _>  HMIN 

p  =  3  corrector  failed  to  converge  for  H  >  HMIN 

p  =  4  corrector  failed  to  converge  for  a  first  order  method 

p  =  5  error  return  from  subroutine  NUITSL 

p  =  6   error  return  from  subroutine  DERVAL 
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MAXDER- input.   Maximum  order  method  which  should  be  used.   The 
highest  order  possible  is  six. 

IPRT-  input.   Print  control  indicator. 
<_  0  ,   no  print  from  LDASUB 

>  0  ,  at  each  step,  print  number  of  steps,  number  of  Jacobian 
evaluations , current  order  being  used,  stepsize  for  next  step, 
current  time,  and  current  values  of  the  dependent  variables. 

H  -  input  and  output.  On  initial  call  it  is  an  estimate  of  the 

timestep.   During  execution  it  is  updated  to  the  current  value, 
and  on  return  contains  the  stepsize  to  be  tried  for  the  next 
step.   The  input  value  need  not  be  accurate.   It  is  better  to 
underestimate  than  to  overestimate  the  initial  value.   The 
stepsize  and  order  are  varied  to  meet  the  local  error  tolerance 
specified.   The  user  does  not  normally  change  the  stepsize 
between  entries  to  SDESOL. 

HMIN-  input.   The  minimum  stepsize  to  be  allowed. 

HMAX-  input.   The  maximum  stepsize  to  be  allowed. 

RMSEPS-input.   The  error  test  constant.   The  values  of  the  relative 
local  errors  must  have  root-mean-square  norm  less  than  RMSEPS. 

W  -  an  array  of  auxiliary  storage  required  by  LDASUB.   This  array 
must  contain  a  total  number  of  locations  equal  to  the  sum  of 
(i)   13*NY  +  5*NL  for  arrays  used  in  LDASUB,  (ii)  storage 
for  the  Jacobian  matrix,  and  (iii)  any  locations  used  in 
processing  the  Jacobian,  e.g.,  scratch  storage  used  by  an 
equation  solver. 
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4.1.2   Subroutine  LDASUB 

This  subroutine  is  the  basic  integration  routine  and  performs  the 
process  in  essentially  the  same  manner  as  subroutine  DFASUB.   A  brief 
description  is  given  in  Section  2  and  differences  between  this  routine 
and  DFASUB  are  outlined  in  Section  3.   Parameters  in  this  routine  in 
which  the  user  may  be  interested  are  stored  in  the  W  array,  an  argument 
of  subroutine  SDESOL. 

YMAX  -  array  of  maximum  magnitudes  of  the  independent  variables, 

y.,   up  to  the  current  time  (or  one,  if  less  than  one). 

This  is  stored  beginning  at  location  7*NY  +  NL  +  1  of 

the  W  array. 
ER    -  This  is  the  array  of  differences  between  the  predicted 

and  corrected  values  of  the  variables,  y.  ,  and  is 

proportional  to  the  estimated  error.   This  array  is 

stored  beginning  in  location   8*NY  +  NL  +  1   of  the  W 

array. 
This  subroutine  incorporates  a  restart  capability.   In  order  to 
restart  from  a  previous  point  without  beginning  again  with  a  first  order 
method,  it  is  necessary  to  have  saved  a  number  of  variables  internal  to 
LDASUB,  and  then  restore  them  before  calling  SDESOL  again.   To  save  the 
internal  parameters,  the  user  calls  subroutine  LDASAV(SAV).   Here  SAV 
is  an  array  of  length  29  in  which  the  values  to  be  saved  will  be  stored. 
In  addition  to  SAV,  the  user  must  also  save  a  number  of  arrays  in  the 
calling  sequence  of  LDASUB,  and  this  is  most  easily  accomplished  by 
saving  the  W  array  in  the  calling  sequence  for  SDESOL.   Once  these 
arrays  have  been  saved,  along  with  the  other  simple  parameters  in  the 
calling  sequence  (  Y   and  YL  need  not  be  saved) ,  the  user  is  free  to 
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use  the  package  to  solve  a  different  problem,  or  to  terminate  the  computer 
run,  to  be  restarted  later. 

At  the  time  the  problem  is  to  be  restarted,  the  user  calls  sub- 
routine LDARST(SAV),  where  SAV  is  the  array  of  values  obtained  previously 
by  calling  LDASAV.   This  restores  internal  values  in  LDASUB.   The  user 
then  calls  SDESOL  with  the  same  simple  parameters  and  the  W  array  as 
before,  except  that  JSKF  =  -  1  and  a  new  end  time,  TEND,  is  provided. 
Restoration  of  values  (including  Y  and  YL  )  in  LDASUB  is  completed 
and  solution  of  the  problem  resumes. 

If  the  user  desires  to  change  the  error  tolerance,  number  of 
variables  in  the  error  test,  or  maximum  order  to  be  used,  the  user 
must  make  a  new  initial  call  to  SDESOL,  that  is,  set  JSKF  =  0  . 

4.1.3  Subroutine  COPYZ 

This  subroutine  simply  transfers  the  contents  of  one  array  into 
another  array. 

4. 2   Supporting  Subroutines 

This  group  of  subroutinesmust ,  at  least  in  part,  be  supplied  by 
the  user.   The  user  must  supply  at  least  one  subroutine,  DIFFUN.   For 
better  efficiency,  the  user  should  supply  a  subroutine,  JACMAT,  to 
explicitly  evaluate  the  Jacobian,  although  a  version  which  approximates 
the  Jacobian  by  numerical  differencing  is  given  in  the  appendix.   To 
take  advantage  of  sparsity  or  other  features  of  his  problem,  the  user 
will  need  to  supply  the  subroutine  NUITSL  to  solve  the  systems  of 
equations  (2.4).   For  certain  problems  the  user  may  have  to  supply 
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subroutine  DERVAL  to  calculate  the  Initial  values  of  the  derivatives. 
We  discuss  the  requirements  of  these  subroutines  in  turn. 

A. 2.1  Subroutine  DIFFUN 

This  subroutine  simply  evaluates  the  equations  (2.1)  at  a  given 
time  and  yalues  of  y  ,  y  ,  and  V   .   Other  parameters  in  the  function 
definition  must  be  transmitted  from  the  calling  program  via  COMMON  or 
some  other  device,  determined  by  the  user. 

The  calling  sequence  is 
CALL  DIFFUN  (Y,  YL,  T,  HINV,  DY)  ,  where  the  parameters  are  defined  as 
follows. 

Y     -  input.   Same  as  in  SDESOL.   This  array  contains  the  current 
values  of  the  variables  y.   and  their  (scaled)  derivatives, 

YL    -   input.   Same  as  in  SDESOL.   This  array  contains  the  current 
values  of  the  linear  variables. 

T     -   input.   Current  time. 

HINV  -   input.   1/h  ,  where  h   is  the  current  stepsize. 

DY    -  output.   Array  of  function  values. 

4.2.2   Subroutine  JACMAT 

This  subroutine  evaluates  the  Jacobian  matrix  J  ,  equation  (2.5) 
at  the  given  time  and  current  values  of  the  dependent  variables,  order, 
and  stepsize.   A  version  of  JACMAT  which  approximates   J  by  numerical 
differencing  is  given  in  the  appendix.   For  maximum  efficiency,  the 
user  should  supply  the  explicit  representation  of  the  Jacobian.   Because 
the  Jacobian  is  used  to  solve  for  the  quasi-Newton  iterates,  it  is  not 
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necessary  for  the  Jacobian  to  be  exact.   Thus  the  user  should  consider 
the  possibility  of  approximations  which  reduce  the  total  number  of 
computations  in  this  step,  with  due  regard  for  the  fact  that  a  smaller 
timestep  may  be  required  to  obtain  convergence  of  the  corrector  within 
three  iterations. 

The  calling  sequence  for  this  subroutine  is 
CALL  JACMAT  (Y,  YL,  T,  HINV,  A2 ,  N,  NY,  EPS,  DY,  Fl ,  PW) ,  where  the 
parameters  are  defined  as  follows. 

Y     -   input.   Same  as  in  SDESOL,   Y  contains  the  current  values 

of  the  variables  y.   and  their  (scaled)  derivatives. 
YL    -  input.   Same  as  in  SDESOL.   This  array  contains  the  current 

values  of  the  linear  variables. 
T     -   input.   Current  time. 

HINV  -  input.   1/h  ,  where  h  is  the  current  stepsize. 
A2    -   input.   The  constant  aQ/$   from  LDASUB. 
N     -   input.   Total  number  of  variables. 

NY    -  input.   Number  of  differential  and  nonlinear  variables. 
EPS   -  input.   Error  constant  from  LDASUB,  vfr  'RMSEPS. 
DY    -  input.   Array  of  current  function  values. 
Fl    -  scratch  array  of  N  locations  available  for  use  by  this 

routine. 
PW    -  output.   The  Jacobian  matrix  J  ,  or  an  approximation, 
calculated  in  JACMAT  and  returned  to  calling  program. 
This  matrix  is  used  in  subroutine  NUITSL  and  the  storage 
mode  must  agree  between  the  two  subroutines. 
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4.2.3  Subroutine  NUITSL 

This  subroutine  solves  the  equations  (2.4)  for  the  quasi-Newton 
iterates.   This  subroutine  will  normally  be  supplied  by  the  user, 
although  versions  which  solve  the  system  by  elimination  methods  and 
iterative  methods,  respectively,  are  given  in  the  examples  in  the 
appendix.   This  subroutine  will  often  be  modified  or  replaced  by  the 
user  to  take  advantage  of  sparsity  or  other  features  of  his  problem  in 
connection  with  JACMAT,  of  course. 

The  calling  sequence  for  this  subroutine  is 
CALL  NUITSL  (PW,  DY,  Fl,  N,  NY.  EPS,  YMAX,  NEWPW,  KRET) ,  where  the 
parameters  are  defined  as  follows. 

PW      -  input.   The  Jacobian  matrix  J  computed  in  JACMAT. 

DY      -  input.   Right  hand  side  of  the  linear  system  to  be 
solved. 

-  output.   The  solution  is  returned  in  the  array  Fl  . 

-  input.   Total  number  of  variables. 

-  input.   Number  of  differential  and  nonlinear  variables. 

-  input.   Error  constant  from  LDASUB,   vM~  •  RMSEPS  . 


Fl 

N 

NY 

EPS 

YMAX 


-  input.   Array  of  maximum  magnitudes  of  y.   up  to  the 


current  time  (or  one  if  maximum  magnitude  is  less  than  one) 
NEWPW   -  input.   Indicates  whether  a  new  J  matrix  has  been 
computed  since  the  last  entry  to  NUITSL. 
-  1,  indicates  this  is  a  new  J  matrix.   If  any 
preprocessing,  such  as  LU  decomposition  is  to  be  done, 
the  preprocessing  should  be  done  and  NEWPW  set  to  zero. 
=  0,  indicates  the  J  matrix  is  the  same  as  on  the  previous 
entry  to  NUITSL. 
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KRET   -  output.   Return  indicator. 
=  0  ,  normal  return 
=  1  ,  error  return,  solution  of  equations  not  obtained. 

Note  that  the  parameters  EPS  and  YMAX  are  useful  if  an  iterative 
method  is  used  for  solutuon  of  the  equations.   Because  the  solution 
represents  corrections  to  the  predicted  value,  and  corrections  to  that, 
the  solution  is  small  compared  to  the  dependent  variable  values.   Hence, 
compared  to  the  YMAX  array,  the  error  tolerance  can  be  fairly  large. 
The  following  convergence  criteria  have  been  used,  with  great  success. 
Let  6u.   denote  the  i —  component  of  the  difference  between  successive 
iterates,  with  u.   being  the  i —  component  of  the  current  iterate. 
Then  the  iteration  is  considered  to  have  converged  whenever 


NY  /  6u.   \2    /      .2 

NY  /    6u.      \2 
±t1\   max(|ui| ,e)  / 


Condition  (i)  requires  convergence  to  2  digits  more  accuracy  than  the 
user  has  asked  for  in  the  solution  of  the  system  (2.1),  relative  to 
YMAX.   Condition  (ii)  requires  the  same  relative  accuracy  in  u.   as  is 
asked  for  by  the  user  in  the  solution  of  the  system  (2.1),  unless  the 
solution  is  smaller  than  e  ,  in  which  case  the  change  is  compared  to 
e  rather  than   lu. I  .   This  avoids  difficulty  if  u.   is  close  to  zero. 
The  e   above  is  EPS  =  vfc-RMSEPS,  where  M  and  RMSEPS   are  inputs  to 
SDESOL.   Two  versions  of  NUITSL  incorporating  iteration  and  this  conver- 
gence test  are  given  in  the  appendix. 
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4.2.4   Subroutine  DERVAL 

This  subroutine  solves  for,  or  otherwise  supplies  the  initial  values 
of  the  derivatives,  and  possibly  other  variables.   In  some  instances 
it  may  need  to  be  supplied  by  the  user.   The  standard  version  of 
DERVAL  given  in  the  appendix  uses  Newton's  method  to  solve  the  first 

m  (=NY)   of  the  equations  (2.1)  for  y(t  )  ,  assuming  values  for  y(t  ) 

3F 
and  V(tn)   have  been  supplied.   To  accomplish  this,  the  matrix  — 

0  3V 

-20  y 
is  needed,  and  this  is  obtained  by  calling  JACMAT  with  h  =  16     , 

A2  =  -1  ,  and  N  =  NY  ,  implying  NL  =  0   for  this  call.   Special  care 

must  be  taken  if  in  fact  NL  is  not  zero  to  assure  that  the  matrix  is 

computed  and  stored  properly.   The  matrix  returned  is  then  16   —  . 

3y 

A  call  to  DIFFUN  yields  the  function  values 

F(y(tQ)  ,  y(tQ)  ,  tQ)  +  PV(tQ)  where  y(tQ)   is  the  current  iterate. 

20 
Multiplication  of  the  function  values  by  16    and  a  call  to  NUITSL 

(again  with  N  =  NY)   gives  the  Newton  iterate.   Of  course,  the  same 

sort  of  special  care  as  necessary  in  JACMAT  is  necessary  in  NUITSL. 

3F 
Obviously  the  above  scheme  cannot  work  if  — -  is  singular,  such 

3y 

as  it  would  be  if  one  of  the  equations  is  algebraic.   In  this  instance 
the  user  must  either  devise  his  own  version  of  DERVAL,  or  supply  the 
values  along  with  a  dummy  version  of  DERVAL.   In  an  extreme  case  the 
user  may  simply  set  initial  derivatives  to  zero.   This  will  provide  a 
poor  predicted  value  on  the  first  step,  and  will  force  an  artificially 
small  timestep  for  the  first  two  steps.   However,  the  overall  penalty 
is  generally  small,  as  appropriate  (corrected)  values  are  computed  at 
the  first  step,  and  after  two  steps  the  program  quickly  increases  the 
timestep. 
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The  calling  sequence  for  this  subroutine  is 
CALL  DERVAL  (Y,  YL,  T,  N,  NY,  DY,  KERET)  ,  where  the  parameters  are 
defined  as  follows. 

Y      -  input  and  output.   Same  as  in  SDESOL.   On  entry  Y(l,i) 

contains  the  initial  values  of  the  variables  y   .   On 

return,  the  values  of  the  derivatives  are  stored  in  Y(2  i)  . 
YL     -  input.   Same  as  in  SDESOL.   This  array  contains  the  initial 

values  of  the  linear  variables. 
T      -  input,   initial  time 
N      -  input.   Total  number  of  variables. 

NY     -  input.  Number  of  differential  and  nonlinear  variables. 
W      -  The  scratch  array  from  SDESOL,  can  be  used  in  any  way 

needed  by  this  subroutine. 
KERET  -  output.   Return  indicater 

=  0  normal  return 

=  1  error  return,  initial  values  were  not  obtained. 
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Appendix  1:   Program  Listings 

The  following  are  listings  of  the  basic  subroutine  package  and 
supporting  subroutines  which  are  of  general  use.   For  simple  problems 
the  user  only  needs  to  supply  a  calling  program  and  a  subroutine,  DIFFUN, 
to  evaluate  the  equations.  Use  of  the  NUITSL  routine  in  computer  facilities 
which  do  not  subscribe  to  the  IMSL  package  will  necessitate  modifications 
to  replace  LUDATF  with  another  LU  decomposition  routine,   and  LUELMF 
with  another  forward  and  backward  substitution  routine. 
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SLBROUTINE  SDESOL  ( Y , YL ,T .TEND, NY  ,NL i M, J3K F  ,MAXDE R  ,  I  PRT  ,H,HMIN,  SOE  ID 

1H*AX,RMSEPS,W)  SCE  20 

C  SCE  30 

c      S0E  t,c 

C  SCE  50 

C      SLBROUTINE  SCESOL  !S  A  DRIVER  ROUTINE  FOR  SUBROUTINE  LCASU6.  SDE  60 

C      ITS  PURPOSE  IS  TO  SET  UP  THE  NECESSARY  REFERFNCE5  TO  A  LARGE  "DE  70 

C      BLOCK  OF  AUMLLARY  STORAGE,  AND   DETAIN  INITIAL  VALUES  OF  SCE  60 

C      DERIVATIVES.  SOE  90 

C      THE  CALLING  SEQUENCE  FOR  SDESOL  IS  SDE  100 

C  SDE  110 

C  CALL  SDESOL(Y,YL,T,TEND,NY,NL,M, JSKF, MAXDER, I PRT, H, HMI N , HMAX, R^SEPS, W ) S CE  120 

C  SDE  130 

C      WHERE  THE  PARAMETERS  ARE  DEFINED  AS  FOLLOW?.  SCE  140 

C  SDE  150 

C  Y        -  ARRAY  DIMENSIONED  <7,NY).   THIS  ARRAY  CONTAINS  T\^=  SCE  160 

C  CEPENDENT  VARIABLES  AND  TH^IR  SCALED  DERIVATIVES.  SDE  170 

C  Y(J+1,I)  CONTAINS  THE  J-TH  DERIVATIVE  OF  THE  I-TH  VARSDE  180 

C  IABLE  TIMES  H**J/ J-FACTORIAl ,  WHERE  H  IS  THE  CURRENT  SDE  190 

C  STEPSIZE.   ON  FIRST  ENTRY  THE  CALLER  SUPPLIES  THE  SDE  20C 

f  INITIAL  VALUES  OF  EACH  VARIABLE  IN  Yd, I).   ON  Sue-  SDE  210 

C  SEQUENT  ENTRIES  IT  IS  ASSUMED  THE  ARRAY  HAS  NOT  SOE  220 

C  BEEN  CHANGED.   TO  INTERPOLATE  TC  NON-MESH  POINTS,  SDE  220 

C  THESE  VALUES  CAN  BE  USED  AS  FOLLOWS.   IF  H  IS  THE  SCE  240 

C  CURRENT  STEPSIZE  AND  VALUES  AT  TIME   T  +  E   ARC  SDE  250 

C  NEEDED,  LET   S  =  E /H   ANC  THEN  SDE  260 

C  SCE  270 

C  JS  SDE  280 

C  I-TH  VARIABLE  AT   T+E  IS   SUM  Y  (  J  +  l ,1 )*S**J  SDE  290 

C  J=0  SCE  300 

C  SDE  210 

C  THE  VALUE  OF  JS  IS  OBTAINED  IN  THE  CALLING  PRCGRAM  SCE  220 

C  BY   JS  =  IAbS( JSKF/10)  SCE  230 

C  YL         APRAY  CF  NL  VARIABLES  WHICH  APPEAL  LINEARLY.  SDE  24C 

C  T        -  CURRENT  VALUE  OF  THE  INCEPENCFNT  VARIABLE  (71*=)  SCE  250 

C  TENC     -  END  TIME  SDE  260 

C  NY       -  NUMBER  OF  DIFFERENTIAL  EQUATICNS  AND  NONLINEAR  SDE  270 

C  VARIABLES.  SDE  280 

C  NL       -  NUMBER  OF  LINEAR  VARIABLES  SDE  29C 

CM-  NUMBER  OF  VARIABLES  INCLUDED  IN  THE  ERROR  TEST  SCE  400 

C  JSKF     -  AN  INDICATOR  JSEC  3CTH  ON  INPUT  AND  OuTPLT  SDE  410 

C  ON  INPUT,  JSKF  =  -1  INDICATES  A  RESTART  CALL  TO  SDE  420 

C  SDESCL.   JSKF  =  0  INDICATES  AN  INITIAL  C»LL  TO  SDE  420 

C  SDESCL.   JSKF  >  0  INDICATES  A  CONTINUATICN  OF  THE  SCE  440 

C  PREVIOUS  CALL  TO  SDESOL.   JSKF  <  -1  MAY  HAVE  RESULTcOSDE  450 

C  FROM  THE  USER  NEGLECTING  TO  TEST  FOR  EPRCR  RETURNS  SCE  460 

C  FROM  SDESCL.   BECAUSE  OF  THIS  POSSIBILITY,  JSKF  <  -1  SCE  470 

C  RESULTS  IN  TERMINATION  CF  THE  RUN  WITH  T\-~  SDE  460 

C  APPROPRIATE  COMMENT.  SCE  490 

C  ON  OUTPUT,  JSKF  CONSISTS  OF  TwO  DIGITS  AND  SIGN,  SCE  500 

C  ♦  OR  -  QP.   Q  IS  THE  ORDER  OF  THE  FORMULA  CURRENTLY  SCE  510 

C  eFING  USED.  P    INDICATES  THE  TYPE  OF  RETURN,  AS  SDE  520 

C  FOLLOWS  SCE  S^C 

C  JSKF  >  6,  P  =  1  IS  THE  NORMAL  RETURN  SDE  540 

C  JSKF  <  0  IS  AN  ERROR  RETURN,  HTH  THE  FOLLOWING  SDE  550 

C  MFAMINGS.  SDE  5fcO 

C  P  =  1   ERROR  TEST  FAILURE  FOR  H  >  HMN  SDE  570 

C  P  =  3   CORRECTOR  FAILED  TC  CONVERGE  FO^  H  >  HMINSDE  580 

C  P  =  4   CORRECTOR  FAILED  TC  CONVERGE  FOR  FIRST  SDE  590 

C  ORDER  METHOD  SDE  600 

C  P  =  5   ERROR  RETURN  FROM  SUBROUTINE  <JU!TSL  SDE  610 

C  P  =  6   ERROR  RETURN  FROM  SUBROUTINE  DE^VAL  SCE  620 

C  MAXDER   -  MAXIMUM  ORDER  DERIVATIVE  THAT  SHOULD  BE  USED  IN  SOE  630 

C  METHOD.   IT  M'JST  BE  NO  GREATER  THAN  SIX.  SCE  640 

C  IPRT     -  INTERNAL  PRTNT  CONTROL  INDTCATCR  FOk  LDASU3.  SCE  65C 

C  IPRT  =  0   NO  PRINT  SDE  660 

C  IPRT  >  0   PRINT  COUNTERS,  STEPSIZE,  CURRENT  TIMESOE  670 

C  AND  VALUES  OF  DEPENDENT  VARIABLES  AT  SDE  680 

C  EACH  STEP.  SDE  69C 

C  H        -  CURRENT  STEPSIZE.   AN  INITIAL  VALUE  ^LST  Br  SUPPLIED  SDE  700 

C  BUT  NEED  NOT  BE  THE  ONE  WHICH  MUST  BE  USED,  SINCE  THESDE  710 

C  SUBRCUTINE  WILL  CHOSE  A  SMALLER  ONE  IF  NECCESSAPY  TO  SDE  720 

C  KEEP  THE  ERROR  PER  STEP  SMALLER  THAN  THE  SPECIFIED  SDE  730 

C  VALUE.   IT  IS  BETTER  TO  UNDERESTIMATE  THE  INITIAL  SOE  740 

C  STEPSIZE  THAN  TO  OVERESTIMATE  IT.   THE  STEPSIZE  IS  SCE  750 
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C                    NORMALLY  NOT  CHANGED  PY  THE  USER.  SDE  763 

C         HMIN     -  MINIMUM  STEPSIZE  ALLOWED  SOE  770 

C         HMAX     -  MAXIMUM  STEP5IZE  ALLOWED  «DE  780 

C         RMSEPS   -  THE  ERROR  TEST  CONSTANT.   THE  CCjCT-mc  AN-SQUARE  OF  SCE  790 

C                   THE  SINGLE  STEP  ERROR  ESTIMATES,   ERU),  QTVIDCC  3Y  SCE  800 

C                    YMAX(I)  =  (MAXIMUM  TO  CURRENT  TIME  DF  Ytl)  )  MUST  3E  SDE  81C 

C                   LESS  THAN  EPS.   THF  STEPSIZE  AND/OR  THE  CRCER  SCE  820 

C                    ARE  VARIED  TO  ACHIEVE  THIS.  SDE  830 

C          W          SCRATCH  STOPAGE  ARRAY.   MUST  EE  AT  LEAST  13*NY  ♦  5*NLSDE  840 

C                   LOCATIONS!  PLUS  THOSE  REQUIRED  FOR  STOPAGE  01=  THE  SDE  850 

C                   MATRIX   PW  (SEE  DESCRIPTION  CF  SUBROUTINE  JACMAT).  SCE  660 

C                    THE  STORAGE  JF  PW  WILL  NORMALLY  REOUIRE  NO  MORE  THAN  SCE  87C 

C                   N**2  +  2*N  LOCATIONS,  AND  IF  C0VPA1T  STORAGc  TECH-  SCE  830 

C                    NIQUES  ARC  USED,  CAN  Be  MUCH  FEWtK.  SDE  890 

C  SCE  900 

c      SDE  910 

DIMENSICN  Y<7,1),  YL(l),  W(l)  SOE  920 

IF  (JSKF.GT.O)  GO  TO  120  SDE  920 

IF  (JSKF.LT.-l)  GC  TO  140  SDE  940 

N  =  NY  +  NL  SCE  950 

IF     (JSKF.LT.O)    GO    TO    110  SCE  960 

C  SDE  97C 

C      IF  THIS  IS  THE  FI«ST  ENTRY,  OBTAIN  VALUES  CF  THE  DERIVATIVE?.  SDE  980 

CALL  DERVAL  ( Y, YL  ,T ,N, NY , W , KRE7R )  SDE  990 

IF  (KRF.TF.NE.O)  GC  TC  130  SDE  1000 

C  SCE  1010 

f      NCW  SET  LP  STORAGE  BLOCKS  IN  THE  W  ARRAY.   THIS  NEEDS  T 0  BE  DONE  SHE  1020 

C      CNLY  INITIALLY  AND  ON  RESTARTS.  SDE  1C20 

C  SDE  10AC 

C                        THE     ARRAY          SAVE          STARTS     AT     LOCATION       1               IN    THE    W    ARPAY  SDE  1050 

C                        THE     ARRAY          YLSV          STARTS     AT    LOCATION       NSVL       I  ISf    THE    W    mRRAY  SCE  1060 

C                       THE    ARRAY         YMAX         STARTS    AT    LOCATION       NrMAX    IN    THE    V,    ARRAY  SCE  1070 

C                        THE     ARRAY          ER                STARTS     AT    LIGATION       NER          IN    THE    W    ARiAY  SCE  1C8C 

C                        THE     ARRAY          ESV            STARTS    AT    LOCATIJM       NESV       IN    THE    W    ARRAY  SCE  1090 

C                        THE    ARRAY          Fl                STARTS    AT    LOCATIJN       NF1          IN    THE    In     AR&AY  SCE  1100 

C                        THE     ARRAY          DY                START?    AT     LOCATION       NCY          IN    THE    W    6R9AY  SCE  1110 

C                       THE    MATRIX       PW               STARTS    AT    LOCATION       NPW         IN    THE    W    ARRAY  SCE  1120 

C  SDE  113C 

110  NSVL  =  7*NY+1  SCE  1140 

NYMAX  =  NSVL+NL  SCE  1150 

NER  =  NYMAX+NY  SDE  1160 

NESV  =  NEP+NY  SDE  1170 

NF1  =  NESV+NY  SDE  1180 

NCY  =  NF1+N  SDE  1190 

NPW  =  NCY+N  SDE  1200 

120  JS  =  JSKF  SDE  1210 

CALL  LDASUB  (  Y,  YL  ,T  ,TFND,  N  ,  vlY,  A,  JS  ,  KF  ,  MAXDEP  ,  I  PRT  ,  t- ,  H^IN,  HM  AX  ,  CCE  1220 

1RMSEPS,W,W(NSVL  )  ,  M  (  NYMAX)  ,  W  (  NER  I  ,  W  (  NE  S  V  I  ,  W  (  N  Fl )  ,  W  (  NC  Y  )  ,  I*  (  MP  W  >  )  SDE  12  20 

C  SCE  1240 

C      CODE  JSKF  ON  RETURN  FROM  LDHSUB  SDE  1250 

C  SOE  1260 

JSKF  ■=  ISIGN(JS*10  +  IABS(KF  ),KF)  SOE  1270 

RETURN  SDE  1280 

130  JSKF  =  -6  SCE  1290 

RETURN  SDE  1200 

140  PRINT  1,  JSKF  >CE  1210 

STOP  iCE  1320 

C  SDE  1220 

C  SCE  1240 

1  FORMAT  COIT  IS  AN  ERROR  TH  ENTER  SDESOL  WITH  JSKF  =  ',110//  SCE  1250 

1  '    RUN  HAS  EEEN  TERMINATED.')  SDE  126C 

END  SDE  1270 


SLBROLTINE  LCASU6  ( Y, YL, T , TEND.M ,NY ,M , J ST  ART ,KFLAG , M A XOR , I PRT f H, 
lHMIN,HMAX,RMSEPS,SAVE,YL3V,YMAX,Ea,ESV,cl,DY,PW) 


LCA  10 

LCA  20 

C                                                                             LCA  20 

C      SLBROUTINE  LCASUb  IS  A  MODIFICATION  OF  SUBROUTINE  DFASUE           LCA  40 

C      kvHICH  IS  DUE  TO  R.  L.  EROWN  AND  C.  W.  GEAR.   DCASUH  IS  DOCUMENTED  LCA  50 

C      IN  THE  REPORT                                                       LDA  60 

C         DOCUMENTATION  FOR  DFASUB —                                       LCA  73 

C         BY  R.  L.  EROWN  AND  C.  W.  GEAR                                    LDA  80 

C         REPORT  LILCDCS-R-73-575.  JUlY  1973                              LDA  90 

C         UNIVERSITY  OF  ILLINOIS  AT  UR EANA -CHAMP  A  I GN                      LDA  100 

C          URBANA,  ILLINOIS   61801                                            LCA  110 
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C  THIS  REPCRT  IS  AVAILABLE  FROM  THE  NATIONAL  TECHNICAL  INF DRM ATI  ON'  LDA  120 

C  SERVICE  CF  THE  'J.  S.  DEPARTMENT  OF  COMMERCE  UNDER  ACCESSION  NUMBERLDA  130 

C  CCO-1469-225 .  LCA  140 

C  LDA  150 

C  THE  MODIFICATION  HERE  IS  DOCUMENTED  IN  THE  REPORT  LDA  160 

C  A  PRCGRAH  FOR  THE  NUMERICAL  SGLUTIOM  OF  LARGd  SPARSE  SYSTEMS  OFLCA  170 

C  ALGEERAIC  AND  IMPLICITLY  DEFINED  STIFF  DIFFERENTIAL  EQUATIJNS  LDA  180 

C  eY  RICFAPC  FRANKE  LCA  190 

C  REPCRT  NPS53FI76051,  MAY  1976  LDA  200 

C  NAVAL  POSTGRADUATE  SCHOOL  LDA  210 

C  MONTEREY,  CALIFORNIA   93940  LCA  220 

C  LDA  230 

C  L0;i  240 

C  LDA  250 

C  ThE  CALLING  SEQUENCE  FDR  LDASUe  IS  LDA  260 

C  .  „„  LDA  270 

C  CALL  LCASUE (Y,YL,T,TEND,LV)4/y,M,JSTAPT,KFLAG,MAX0Rt IPRT,H,HMIN,  LDA  280 

C  HMAX,RMSEPS,SAVE, YLSV,  YMAX",  E"  ,  ESV  ,  F  1 ,  DY,  PW  )  LDA  290 

C  LCA  300 

C  WHERE  THE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS.  LDA  310 

C  Y        -  ARRAY  DIMENSIONED  <7TNY).   THIS  ARRAY  CGNTAINS  THE  LDA  320 

C  DEPENDENT  VARIABLES  AMD  THEIR  SCALED  DERIVATIVES.  LCA  330 

C  Y(J+1,I)  CONTAINS  THE  J-fw  DERIVATIVE  OF  THE  I-TH  VARLDA  340 

C  IABlE  TIMES  H**J/J-FACTORIAL,  aHERE  H  IS  THE  CURRENT  LOA  350 

C  STEPSIZE.   ?N  FIRST  ENTRY  THE  CALLER  SUPPLIES  THE  LCA  360 

C  INITIAL  VALUES  OF  EACH  VARIABLE  IN  Yd, I)  AsD  AN  LCA  370 

C  ESTIMATE  OF  THE  INITIAL  VALUES  CF  THE  CERIVATIVFS  LCA  380 

C  IN  Y(2,I>.   ON  SUBSEQUENT  ENTRIES  IT  IS  ASSUMED  THAT  LCA  390 

C  THE  ARRAY  HAS  NOT  BEEN  CHANGED.   TO  INTERPOLATE  Tn  LCA  400 

C  NON-MESH  POINTS,  THESE  VALUES  CAN  BE  USEC  AS  FOLLOWS. LCA  410 

C  IF  H  IS  THE  CJRRENT  STEPSIZE  AND  VALUE?  AT  TIME  T+E  LDA  420 

C  NEEDED,  LET  S  =  E /H   AMD  THEN  LCA  430 

C  LDA  440 

C  NO  LDA  450 

C  I-TH  VARIABLE  AT   T+E  IS   SUM  Y < J+l  ,  I  )*S** J  LCA  460 

C  J=0  LDA  470 

C  LDA  480 

C  THE  VALUE  OF  NQ  IS  OBTAINED  IN  THE  CALLING  PROGRAM  LDA  490 

C  BY   NQ  =  JSTART.  LDA  500 

C  LCA  51.1 

C  YL          ARRAY  OF  NL  =  N  -  NY  VARIABLES  WHICH  APPEAR  LINEARLY. LCA  520 

C  THE  USER  SU°PLIES  INITIAL  VALUES  FOR  THESE  VAR I ABLc S .L DA  £30 

C  T        -  CURRENT  VALUE  OF  THE  INDEPENDENT  VARIABLE  (TIME)  LOA  540 

C  TEND     -  END  TIME  LDA  550 

C  N          TOTAL  NUMBER  OF  VARIABLES  LDA  560 

C  NY       -  NUMBER  OF  DIFFERENTIAL  EQUATIONS  AND  NONLINFAR  LCA  570 

C  VARIABLES.  LDA  580 

CM-  NUMBER  OF  VARIABLES  INCLUDED  IN  THE  ERROR  TEST.  LDA  590 

C  THIS  NUMBER  CAN  BE  NO  GREATER  THAN  NY.   IF  IT  IS  LDA  600 

C  GREATER  THAN  NV,  NY  VARIABLES  ARE  USED  IN  tj-E  ERROR  LOA  610 

C  TFST.  LCA  620 

C  JSTART   -  INPUT  AND  OUTPUT  INDICATOR.  LCA  630 

C  ON  INPUT  JSTART  HAS  THE  FOLLOWING  MEANINGS.  LCA  640 

C  <0   THIS  INDICATES  A  RE-START  FROM  A  PREVIOUS  LCA  650 

C  POINT  FOLLOWING  TcRMIINATIGN  OF  THF  SUN  OR  LDA  660 

C  SOLUTION  OF  ANOTHER  PRCBLFM  DURING  1  HE  SAMF  LCA  670 

C  RUN.   PARAMETERS  IN  THE  CALLING  SEQUENCE  LCA  680 

C  MUST  HAVC  BEEN  PRESERVED  FROM  THE  PPFVIOUS  LCA  690 

C  USE,  PARTICULARLY  THE  ARRAYS  LCA  700 

C  SAVE,  YLSV,  ESV,  AND  PW.  LDA  710 

C  THESE  ARRAYS  MUST  BE  SAVED  AFTER  A  CALL  LDA  720 

C  TO  SUBROUTINE  LDASAV,  WHICH  ALSC  SAVES  LLA  730 

C  NECESSARY  PARAMETERS  INTERNAL  TC  LDASUB.  LCA  740 

C  =0   INDICATES  AN  INI7IAL  CALL  TO  LDASUB.   THE  LCA  750 

C  ROUTINE  INITIALISES  ITSELF,  SCALES  THF  LDA  760 

C  DERIVATIVES  IN  Y(2,I)  AND  THEN  PERFORMS  THE  LDA  770 

C  INTEGRATION  UNTIL  T  >  TEND.  LCA  780 

C  >0   INDICATES  THE  SOLUTION  IS  TO  BE  CONTINUED.  LOA  790 

C  AFTER  THE  INITIAL  ENTRY  IT  Io  NEITHER  LDA  800 

C  DESIRABLE  NOR  NCCESSAHY  TO  RE-ENTER  WITH  LCA  810 

C  JSTART  =  0,  SINCfc  THIS  RE- INI TI AL I ZES  LOA  82C 

C  THE  CODE,  BEGINNING  WITH  A  FIRST  CRDER  LCA  830 

C  METHOD  AGAIN.  LDA  840 

C  ON  OUTPUT,  JSTART  IS  SET  TO  THE  VALUE  OF  NC,  THE  LDA  850 

C  ORDER  OF  THF  FORMULA  CURRENTLY  BEING  USEC.  LCA  860 


24 


c 

KFLAG    -  THE  COMPLETION  CODE  INDICATE,  Wll   THE  FOLLOWING 

IDA 

570 

c 

MEANINGS 

i.  :a 

680 

c 

+1   THE  INTEGRATION  WAS  SUCCESSFUL 

l  CA 

590 

c 

-1   cRRGR  TEST  EATLUGc  fCR  H  >    HMN 

LCA 

900 

c 

-3   CORRECTOR  F»-lLED  TO  CCNVEnGE  FOR  H  >  HMIfv 

LCA 

913 

c 

-4   CORRECTOR  FAIlEC  TO    CCNVFkGE  FO"  FIRST 

lCA 

920 

c 

ORDER  METHOD 

LCA 

920 

c 

-5   ERROR  RETURN  FRQM  SUBROUTINE  NUITSL 

LCA 

940 

c 

MAX"1 

5    -  MAXIMUM  ORDER  DERIVATIVE  THA~  CHCUuD  OE  LSfcD  IN  THE 

LCA 

9  50 

c 

METHOD.   IT  MUST  BF  NO  GREATER  THAN  SiX.   IF  IT  T c 

I.  CA 

960 

c 

GREATER  THAN  SIX,  THE  MAXIMUM  CRDEP  U'ED  WILL  BE  SIX 

.LCA 

970 

c 

IPRT 

-  INTERNAL  PRTNT  CONTROL  INHCATCrt 

LCA 

980 

c 

-    0   NO  PRINT 

LCA 

9  90 

c 

>  0   PRINT  COUNTERS,  "TEPSIZE,  CURRENT  KME 
AND  VALUES  OF  DEPENDENT  VARIABLES  AT 

LDA 

1000 

c 

LCA 

1010 

c 

EACH  STEP. 

1.0  A 

1020 

c 

H 

-  CURRENT  STEPSIZE.   AN  INlTiAL  VALUE  MUST  BE  SUPPLIED 

LDA 

103J 

c 

BUT  NEED  NOT  BE  THE  1\E    wH  CH  wILL  6 F  USED,  SINCE  THfcLCA 
SUBROUTINE  WILL  CHCCSE  A  SMALLER  ONE  '■  F    NECESSARY  TCL  CA 

1040 

c 

1C5C 

c 

KEEP  THE  ERROR  o£R  STEP  SMALLER  ThAN  THE  SPFCIFIEC 

LCA 

1060 

c 

VALUE.   IT  TS  BCTTER  TO  UNDERESTIMATE  THE  INITIAL 

LCA 

1070 

c 

STEPSIZF  THAN  TC  OVFKEiTIMATE  IT.   THE  STEPSIZF  IS 

LCA 

1080 

r 

NORMALLY  NOT  CHANGED  BY  THE  USER. 

LCA 

1090 

c 

HMINI 

-  MINIMUM  STEPSIZE  ALLOWED 

LCA 

1100 

c 

HMAX 

-  MAXIMUM  STEPSIZE  ALLOWED 

LDA 

1110 

c 

RMSE 

FS   -  THE  ERROR  TEST  CONSTANT.   THE  POOT-ME AN-SQU^E  OF 

LCA 

1120 

c 

THE  SINGLE  STEP  ERROR  ESTIMATES,   ER (  I  )  ,  DIVIDED  BY 
YMAXm  =  (MAXIMUM  TO  CURRCNT  TIME  OF  Y(IM  MUST  PE 

LDA 

1130 

c 

LCA 

1140 

c 

LESS  THAN  RMSEPS.   THE  STEPSTZF  AND/OR  CRDEF  ARE 

LCA 

1150 

c 

VARIED  ~1    ACHIEVE  THIS. 

LDA 

1  160 

c 

SAVE 

-  AN  ARRAY  Oc  LENGTH  AT  LEAST  7*NY 

LCA 

1170 

c 

YLSV 

-  AN  ARRAY  CF  LENGTH  AT  LEAST  NL 

LCA 

1180 

c 

YMAX 

-  A  VECTOR  OF  LENGTH  NY  WHICH  CONTAINS  TEE  MAXIMUM 

LDA 

1190 

r 

OF  EACH  Y  SE^N  SO  FAR.   ON  THE  FIRST  CALL,  THESF  WiLLLCA 

1200 

c 

BE  INITIALIZED  AS  YMAX(I)  =  M t X ( 1 ,  | Y (  1 ,  I  )  I  ) 

LCA 

1210 

c 

ER 

-  A  VECTOR  OF  LENGTH  MY 

LCA 

1220 

cc 

ESV 

-  A  VECTOR  GE  LENGTH  NY 

LDA 

12;0 

^1 

-  A  VECTOR  V  LENGTH  N  =  NY  +  *'L 

LDA 

1240 

c 

DY 

-  A  VECTOR  i"-    LENGTH  N  =  NY  +  NL 

LCA 

1250 

c 

PW 

-  AN  ARRAY  IN  WHICH  THE   J   MA"PIX  CCMt>uTEC 

LOA 

1260 

c 

IN  SUBROUTINE  JACMAT  WILL  BE  STORED.   SIZE  WHICH 

LCA 

12  70 

c 

MUST  9E  ALLOWED  IS  DETERMINED  BY  THF  STCRAGE  TECH- 

LCA 

1280 

c 

NIQUE  USED  FOR  IT,  BUT  NORMALLY  WON'T  BE  MORE  THAN 

LCA 

129C 

c 

U**l    +  2*N  LOCATIONS,  THE  LAT"ER  2*N  BEING  REQUIRED 

LCA 

1200 

c 

BY  THE  LINEAR  EQUATION  SOLVER. 

LCA 

1210 

c 

LOA 

12  20 

c 

-LCA 

1330 

OIMENSI 

IN  Y(7,l),  YL(1),  SAVE(7,l),  YMAX(l),  ER(1),  YLSV<1),  P 1 ( 1 ) t 

1340 

1.  PERK 

6,2),  CUF(21),  ESV(l),  DY<1),  PHI),  SAV(l).  A(29) 

ENCE  (A(8),BND),  (  A  (  9  ) ,  BR  )  ,  (A(10),F>,  (  A  (  1 1 )  .  EDW  \i )  , 

=  NQ1).  ( A(  13) ,FNQ2 ) ,  < A { 14 ) , ENP3 ) ,  ( A ( 15 ) , EPS )  ,  (A(i6),EUP 

rr>N«:K-)|  < A ( 18 ) , PEP SH  )  ,  (  A  (  19  )  ,  I  D1U6  )  ,  ( A  (  20  )  ,  I  kEVAL  )  , 

K),  (A(22) ,LCOPYL) .  ( A ( 2 3  )  ,  LCOP Y Y  )  ,  I A ( 24)  , MA \rc R ) , 

Ml),  (A(26),NL),  <A(27),NQ),  (A(28),KS>,  (A(29),Nw) 

LDA 

1250 

ECUIV/SL 

LCA 

12  60 

1<A(12), 

)LDA 

13  70 

2,(A( 17) 

LC« 

1380 

3  {/l  <21 )  . 

LCA 

1390 

4( A(25), 

LCA 

1400 

c 

LCA 

1410 

c 
c 

1420 
14  2C 

LDA 

c 

TEE  COE 

FFICIENTS  IN  THE  PERT  ARRAY  ARE  USED  FOR  tRPCR  TESTING  AND 

LCA 

1440 

c 

CHANGING 

LDA 

14  50 

c 

LCA 

1460 

c 

1470 
1480 

CATA  PE 

RT/4.,9.,16.t25.,36.,49.,9.,16.,25.,:-c.,49.,64.,l.,l.,.25, 

LQA 

12.7889E 

-2,1.7C56  9E-3,6.8  3929E-5/ 

LCA 

1490 

c 
c 

1500 
1510 

LDA 

c 

THE  ENTR: 

LCA 

1520 

c 

STABLE 

METHODS  USED  IN  THIS  PROGRAM  AND  ARE  TO  3c  THE  MACHINE 

LCA 

1530 

c 

PBECI3I 

ON  ECUIVALENTS  OF  THE  FOLLOWING  CONSTANTS. 

LCA 

1540 

c 

LOA 

1550 

c 

-1 

LCA 

1560 

c 

-2/2  , 

-1/2 

LCA 

1570 

c 

-11/6  , 

-1  ,  -1/6 

LDA 

1580 

c 

-25/12 

,  -35/24  ,  -5/12  ,  -1/24 

LCA 

1590 

c 

-137/60 

,  -15/8  ,  -17/2t  ,  -1/8  ,  -1/120 

LCA 

16C3 

c 

-147/60 

,  -2C3/90  ,  -49/48  ,  -35/144  ,  -7/240  ,  -1/720 

LCA 

1610 

25 


C  LDA  1620 

c  LDA  lfc30 

DATA    CO F/-1..-1. 5, -.5, -I. 833333,-1. ,-.1666 667, -2. 083 3 2 3, -1.453333, LCA  1640 

1-.4 1666 6  7, -.041 b6 6 67, -2. 2 83 3 33, -1.8 75 ,-. 7082333, -. 125  •-. 003 333233* LCA  16  50 

2- 2.  45, -2.  2  5  55  56, -1.0  206  33, -.2430556, -.02  9  1666  7, -.00 13  8a  8  89/  LD.t  16  60 

IF     (JSTART)     100,110,150  LCA  1670 

c  Lr,A  lfc80 

C  IF    THIS     15    A    RESTART    ENTRY,     RESTORE    Y    A^D    YL    FROM     IHE    SAVF     ANC  LCA  1690 

C  YLSV    ARRAYS,    KHEP  E    THEY    l*  =  RE    SAVED    BY    A    PREVIOUS     CHI    TC    LDASAV.       LDA  1700 

C  LDA  1?10 

100    CALL    C0PY2     (Y,3AVE,LC0PYY)  LCA  1720 

CALL    COPYZ     (YL, YLSV, LCOPYL  )  LCA  172C 

GC    T0    15C  LDA  1740 

r  LDA  1750 

C  IF    THIS    IS    THE    FIRST    CALL,     INITIALIZE    Y.XAX,     SCALE    DERIVATIVES,    ANDLDA  1760 

C  INITIALIZE    INDICATORS    AND    SET    ORDER    TO    ONE.  LDA  1770 

C  FGR    DOUeLE    PRECISION,    SET    LCOPYL    =    14*NY    AND    LC3PYL    =    2*NL     IF  LD*  1780 

C  SLBROUTINE    CCPYZ     IS    IN    SINGLE    PRECISION.  LDA  179C 

c  LDA  180o 

110    NL    =    N-NY  LDA  161C 

LCOPYY    =    7*NY  LCA  1820 

LCOPYL     »    NL  LCA  1830 

vi  =  mincjcnyj  lda  1840 

EPS  =  SCRT(FLOAT(M)  )*RMSEPS  LDA  1850 

MAXCER  =  VIN0(MAXCR,6)  LCA  I860 

IF  (IPPT.LE.C)  GO  TO  120  LDA  1570 

PRINT  3,  N,ia,RMSEDS,TEND,H  LCA  1880 

PRINT  4  LDA  1890 

120  NS  -    0  LDA  1900 

NW  =  0  LDA  1910 

C  LCA  1920 

CO  130  J=1,NY  LDA  1930 

YfAX(J)  --     &l*AXl(l.,AbS(Y<  1,  J)  )  )  LDA  1940 

130  Y(2,J)  =  Y(2,J)*H  LDA  1950 

C  LDA  1960 

NC  =  1  LDA  1970 

eR  =  1.  LDA  1980 

ASSIGN  190  TC  IRET  LCA  1990 

c      L0A  2000 

f  SET    COEFFICIENTS    FOR    THE    ORDER    CURRENTLY    BEING    USED.  LDA  2010 

C  E    IS    A     TEST    FOR    ERRORS    OF    THE    CURRENT    ORDER    NO  LDA  2020 

C  EUP    IS    TC    TEST    FOR    INCREASING    THE    ORDER,     ECWN    FOR     DECREASING    THE       LDA  2030 

C  CPCER.  LCA  2040 

c  LDA  2050 

140    K    =    NC*(NQ-l>/2  LDA  2060 

CALL    COPYZ     (A(2),C0F(K+1)  ,N0>  LDA  2070 

K    =    NG+1  LCA  2080 

ICOUB    =    NC  LCA  2090 

ENC1    =    .5/NC  LCA  2100 

ENC2    =    .5/K  LCA  2110 

ENC3    =     .5/<NC+2)  LCA  2120 

PEPSH    =    EPS**2  LDA  2130 

E    =    PFRT(KC,1 l*PFPSH  LDA  2140 

cLP    =    DERT(NC,2)*PEPSH  LCA  2150 

EQWN    =     PERT(NC,3)*PEPSH  LDA  2160 

BNC    =     <EPS*E,NQ3)**2  LCA  ' 

UEVAL    =    1  LDA 


UEVAL    =    1 

GC    TO    IRET,     (190,200,490,570) 

IF     (H.FC.HNEW)     GO    TO     190 


2170 

2180 

LCA    2190 

150    IF     (H.FC.HNEW)    GO    TO    190  LCA    2200 

-LCA    2210 

C               IF    CALLER    HAS    CHANGED    H,     RESCALE    DERIVATIVES    TO    REFLECT    THAT    FNFW    LCA    2220 
C               VsAS    USEC    CN    THE    LAST    CALL.  LCA    2230 

c  LDA    2  240 

R  =  H/HNEH  LDA  2250 

ASSIGN  190  TO  IRET  LCA  2260 

GC  TO  61C  LCA  2270 

c      LDA  2  280 

C      SET  JSTART  TC  NQ,  THE  CURRENT  ORDER  OF  THE  FETH°D,  B=F3P=  EXIT,    LDA  2290 
C      AND  SAVE  THE  CURRENT  STEPSIZE  IN  HN^W.  LCA  2200 

C      LCA  2310 

160  JSTART  =  NQ  LDA  2220 

FNEW  =  H  LCA  2330 

RETURN  LCA  2340 

170  NS  =  NS+1  LCA  2250 

IF  (IPRT.LE.O)  GO  TO  180                                             LCA  2360 
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LCfi  2270 

PRINT  OAT*  IP  DESIRED  BY  USES  LDA  2280 
Lr,A  2390 

PRINT  1,  NS ,NW,NQ,H,T,(Y< 1,1 ).I=1,NY)  LCA  2400 
IF  (NL.GT.O)  PRINT  2,  < YL ( I  )  , I -1, NL )  LDA  2410 
180  CONTINUE  LCA  2420 
IF  (KFLAC.LT.O)  GO  TO  160  LCA  242C 
IF  (T.Gfc.TcNC)  GO  TO  160  LDA  2440 
LOtf  24  50 

TAKE  ANCTHtR  STEP  IF  T  <  TEND  LCA  246G 
L  |JA  2 4  70 

JSTART  =  1  LCA  2480 
LDA  2490 

SAVE  DATA  FCP  TRTAL  WITH  A  SMALLER  TIMESTEP  IF  THIS  STEP  FAILS  LCA  2500 
LCA  2  510 

190    CALL    COPYZ     (SAVE, Y,LCOPYY)                                                                                                             LCA  2520 

CALL    COPYZ     (YLSV,YL,LCOPYL)                                                                                                          LCA  2520 

RACUM    =     1.                                                                                                                                                               LCA  2540 

KFLAG     =    1                                                                                                                                                                  LDA  2550 

HCLD    =    b                                                                                                                                                                     LCA  2560 

NCCLD    =    NQ                                                                                                                                                       LCA  2570 

TCLD    =    T                                                                                                                                                            LCA  2580 

2CC    T    =    T*»                                                                                                                                                               LCA  2590 

UNV    =     l./H                                                                                                                                                    LCA  2600 

c               LCA  -tl0 

C                CCMPUTE     FRFClCfED    VALUES     BY     5FFECTTVELY    MULTIPLYING    DcRlV^TIVE             LCA  2620 

C                VECTCR     eY     PASCAL    TRIANGLE     ^A~PIX                                                                                                 LDA  2620 

C                LCfi  2  6  40 

C                                                                                                                                                                                                          LCA  2650 

DC     210    J=2,K                                                                                                                                                          LCA  2660 

J2     =    K+J-l                                                                                                                                                               LCA  2670 

C                                                                                                                                                                                               LDA  2680 

DO     210     J1=J,K                                                                                                                                                       LCA  2690 

J2    =    J3-J1                                                                                                                                                       LCA  2700 

C                                                                                                                                                                                               LCA  2710 

DC  210  1  =  1, NY                                                       LCA  2720 

210  Y4J2,I)  =  Y(J2,I )*Y( J2+1, I )                                         LCA  2730 

C                                                                         LOA  2740 

C                                                                         LCA  2750 

DC  220  1  =  1, NY                                                          LCA  2760 

220  ER(I  )  =  C.                                                          LCA  2770 

LCA  27SO 
LDt  2  790 

DC  U°  TC  THREE  COO  ECTOR  ITERATIONS.  CONVERGENCE  IS  :6TA[NEC  WHENLO*  2800 
CHANGES  ARE  LESS  THAN  6ND  WHICH  IS  DEPENDENT  DN  THE  EP-SCR  TEST  LCA  2810 
CCNSTANT.  ThE  SUM  OF  CORRECTIONS  IS  ACCUMuL^ED  IN  Efi(I).  IT  IS  LCA  2820 
ECUAL  TO  ThE  K-TH  DERIVATIVE  GF  Y  TIMES  H**K / ( K-F AC r0 R I AL * A ( K )  ) ,  LC*  2830 
AND  THUS  IS  PROPORTIONS  TO  THc  ACTUAL  ERRORS  TO  THE  L3i»=ST  PO*FR  LCA  2640 
OF  H  PRESENT,  WHICH  IS  H**K.  LCA  2850 
LDA  2  860 

LDA  2670 
DC  270  L=l,2  LCA  2680 
CAlL  DJPFUN  (Y, YL  ,T,HINV, OY)  LCA  2890 
IF  (IUFVAL.LT.1)  CO  TO  230  LDA  2903 
L0A  2  91C 

IF  THERE  HAS  BFE\  A  CHANGE  CF  3RDER  30.  THERE  HAS  3EEN  TROUPLF  LCA  2920 
WITH  CONVERGENCE,  PW  IS  R E -t V ALUATED  PRIOR  TO  STARTING  TH-  LCA  2930 
CCRPECTCR  ITERATICN.  IWEVAL  IS  THEN  SET  TO  -1  AS  AN  INDICATOR  LCA  2940 
THAT  IT  HAS  EEEN  DONE.  NEWPw  IS  SET  ^jNZERC  TO  i^iiiCt'z  TO  lCA  2950 
cLBROUTINE  NUITSL  THAT  A  NEW  PW  HAS  3fcrN  F0~v!Dc3.  LCA  2960 
LCA  2  9  70 

CALL  JACNAT  ( Y , Yl , T , HI NV , A ( 2 ) , N , NY , E PS , DY , F  I . PW )                    LCA  2980 

KFLAG  =  1                                                              '-CA  2990 

IWEVAL  =  -1                                                         LCA  3C00 

NW  =  NW+1                                                              LCA  2010 

NEWPW  =  1                                                              LOA  3020 

230  CALL  NUITSL  ( PW , D Y  ,  F 1 , N ,MY , E PS , YM AX , NEWP W , KRH ct )                    LCA  3020 

IF  (K3R6T.NE.0)  GC  TO  600                                             LCA  3040 

IF  (NL.LE.O)  GO  TO  250                                              LCA  3050 

LOA  3060 

CO  240  I"U,NL                                                       LCA  3070 

240  YL(I)  =  YL (I )-Fl< I+NY)                                                LOA  3080 

LDA  3090 

25C  CCNTINUE                                                            LCA  3100 

3EL  *  0.                                                            LOA  2110 
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C  LCA  3120 

DC  260  1=1, NY  LOA  3130 

Yd,  I  )  =  Y(1,I)-F1(I)  LD6  3140 

Y(2,I)  =  Y(2,I)*A(2)*Fll I )  IDA  3150 

EP (I )  =  ER<!)+F1( I)  LCA  3160 

DEL  =  DEL-MF1  (I )/ AM AX1 ( YMAX ( I  ),ABS( Y( 1, I ) ) )  )**2  LCA  3170 

260  CONTINUE  LDA  3180 

C  LCA  3190 

IF  (L.GE.2)  ER  =  AMAX1 ( .9* BR , DEL/DEL  1  )  LDA  3200 

DELI  =  CEL  LDA  3210 

IF  <AMIN1(DEL,BR*DEL*2.).LE.BND1  GO  TO  330  LDA  3220 

270  CONTINUE  LDA  3230 

C  LDA  3240 

C      LD£  3  250 

C      THE  CCRRECTIOR  ITERATION  FAILED  TO  CONVERGE  IN  3  TRIES.   VARIOUS  LCA  3260 

C      POSSIBILITIES  ARE  CHECKED  FOR.   IF  H  IS  ALREADY  HMIN  AND  PW  HAS  LCA  3270 

C      ALREADY  BEEN  RE-EVALUATED.  A  NO  CONVERGENCE  ;XIT  IS  TAKEN.  LDA  3283 
C      OTHERWISE  THE  MATRIX  PW  IS  RE-EVALUATED  AND/OR  (IN  THAT  ORDER)  THELDA  3290 

C      STEP  IS  REDLCED  TO  TRY  AND  GET  CONVERGENCE.  LDA  3300 


T  =  TOLC 

IF  (IWEVAL)  280,300,290 
280  IF  (H.LE.HMIN*1. 00001)  GO  TO  310 

290  RACUM  =  PACLM*.25 
300  CONTINUE 

GC  TC  560 
310  KFLAG  =  -3 


C 


•LOA  3310 


C      RESTORE  Y  ANC  YL  AFTER  CONVERGENCE  FAILURE 
C      

320  CALL  COPYZ  (  Y  ,S A VE , LCOPYY ) 
CALL  COPYZ  (YL,YLSV,LCOPYL) 
H  =  HOLC 
NC  =  NOCLO 
GO  TO  170 


C      THE  CORRECTCP  CONVERGED,  SO  NOW  THE  ERROR  TEST  IS  MADE. 
C      

330  C  =  0. 


DC  340  1  =  1, M 

YN  =  AMAX1(AES(Y( 1,  I) ),YMAX( I ) ) 

340  D  =  D*< EP( I )/YM)**2 

UEVAL  =  0 

IF  (C.GT.E)  GO  TO  380 


C  THE  ERRCP  TEST  IS  OKAY.  SO  THE  STEP  IS  ACCEPTED.   IF  IQOUfi 

C  NCW  BECOMES  NEGATIVE,  A  TEST  IS  MADE  TO  SEE  IF  THE  STEP  SIZE 

C  CAN  BE  INCREASED  AT  THIS  ORDER  OR  ONE  HIGHER  OR  ONE  LCWER. 

C  THE  CHANGE  IS  MADE  ONLY  IF  THE  STEP  CAN  BE  INCREASED  PY  AT 

C  LEAST  101.   ICOUB  IS  SET  TO  IMC  TO  PREVENT  FLRTHER  TESTING 

C  FOR  A  WHILE.   IF  NO  CHANGE  IS  *ADE,  IDOUB  IS  SET  TO  9. 


IF  (K.LT.3)  GO  TO  360 

CC  350  J=3,K 

DO  3  50  1=1, NY 
350  Y(J,I)  =  Y(J,I)+A(J)*ER(I) 

360  KFLAG  =  1 

ICOUB  =  IDCLE-1 

IF  (IDOLB)  410,370,510 

370  CALL  COPYZ  (ESV,ER,M1) 
GC  TO  510 


C  THE  FRRCP  TEST  FAILED.   IF  JSTART  =  0,  THE  DERIVATIVES  IN  THE 

C  SAVE  ARRAY  ARE  UPDATED.   TESTS  ARE  THEN  MACE  T0  FIX  THE  STEPSIZt 

C  AND  PERHAPS  RECUCE  THE  ORDER.   AFTER  RESTORING  AND  SCALING  THE 

C  Y  VARIABLES,  THE  STEP  IS  RETRIED. 


380  IF  (JSTART. GT.O)  GO  TO  400 
DO  390  1=1, NY 


LDA 

3320 

IDA 

3  3  30 

LOA 

3340 

LDA 

3  3  50 

LCA 

3360 

LDA 

3370 

LDA 

3380 

LDA 

3390 

LDA 

3400 

LDA 

3410 

LDA 

3420 

LCA 

3430 

LDA 

3440 

1  OA 

3  4  50 

LCA 

3460 

LDA 

3470 

LDA 

3480 

LCA 

3490 

LDA 

3  500 

LDA 

3  510 

LCA 

3520 

LDA 

3  5  30 

LDA 

3  540 

LDA 

3550 

LDA 

3560 

LDA 

3570 

LDA 

3580 

LDA 

3590 

LDA 

3600 

LDA 

3610 

LDA 

3620 

LCA 

3630 

LDA 

3640 

LDA 

3650 

LDA 

366C 

LDA 

3670 

LDA 

3680 

LDA 

3690 

LDA 

3700 

LDA 

3710 

LDA 

3720 

LDA 

3730 

LDA 

3740 

LCA 

3750 

LDA 

3760 

LDA 

3770 

LDA 

3780 

LDA 

3790 

LDA 

3800 

LDA 

3810 

LDA 

3820 

LDA 

3830 

LDA 

3840 

LDA 

3850 

LCA 

3860 
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390  SAVE(2,I)  =  ><2,I  I 


400 
410 

420 

4  30 

440 

450 

460 

470 


KFLAG  = 

IF  (H.L 

T  =  TOL 

IF  (KFL 

PP2  =  ( 

L  =  0 

IF  (NG.LE.l »  GO  TC  4  30 

C  =  0. 


KFLAG-2 
E.HMIN)  GO  TO  5  50 
C 

AC.LE.-5)  GO  TO  530 
C/E)**ENQ2*1.2 


CC  420  J=i,M 

YM  =  AMAXl(AeS(Y( 1,J)  ),YMAX( J)  ) 

C  =  0-M  V(K,  J>/YMJ**2 


PRl  =  ( 

IF  (PRl, 

PR2  =  P! 

L  =  -I 

IF  (KFLAG.  L7.0. OR. NQ.GE.MAXD6R )  GO  TO  450 

C  =  0 


C/EDWN)**ENQ1*1.3 
..GE.PP2)  GO  TO  430 
•RI 


DO  440 
Yf  =  AM 
D  =  D+< 

PRl  =  ( 
IF  (PRl 
PP2  =  P 
I  =  1 
ft  =  1./ 
IF  (KFL 
ICOUB  = 
GO  TO  5 
NEWC  = 
K  =  NEW 
IF  (NEW 
Rl  =  A( 

CC  470 
Y(K, J) 


J«1»FX 

AX1 (AES(Y( 1,J)  ),YMAX<  J)] 

(ER(JI-ESV< J) >/YM)**2 


C/E 
.GE 

Rl 

ANA 
AG. 

9 
10 
NC  + 
C-H 
C.L 
NEW 

J  =  l 
=  E 


UF )**E^G3+1.4 
.FR2)  GO  TO  450 


Xl(PR2t  1.5-5) 

LT.O.OR.R.GE.i.l)  C-0  TO  460 


E.NQ)  GO  TO  480 
Q)/FLOAT(NEWQ) 

.AY 

R(J)*R1 


430  CCNTINUE 


490 
50C 


510 
520 


IF  THF  STEP  WAS  OKAY,  SCALE  THE  Y  VARIABLES  IN  ACCORDANCE 
WITH  ThE  NEW  VALUE  OF  H.   IF  KFLAG  <  0,  HOWEVER,  USE  THE 
SAVED  VALUES  (IN  SAVE  AND  YLSV).   IN  EITHER  CASE,  IF  THE  ORJER 
HAS  CHANGED  IT  IS  NECESSARY  TO  FIX  CERTAIN  PARAMETERS  BY  CALLING 
THE  PROGRAM  SEGMENT  AT  STATEMENT  NUMBER  140. 

ICGUB  =  NC 

IF  (NEwC.EQ.NC)  GC  TO  490 

NC  =  NEWC 

ASSIGN  4S0  TC  IRET 

GO  T"  140 

IF  (KFLAG. GT.O)  GO  TO  500 

RACUM  =  PACUN*R 

GC  TO  560 

R  =  AMAX1(AMIN1(HMAX/H,P) ,HMIN/H) 

F  =  H*R 

IfcEVAL  =  1 

ASSIGN  510  TC  IRET 

GC  TO  610 

CO  520  1=1, H 

YMX(I)  =  AKAX1(AES(Y(  1,1  ))  ,  YMAXU  )  ) 


GC  TO  170 


THE  ERROR  TEST  HAS  NOW  FAILED  THREE  TIMES,  SC  THE  DERIVATIVES  ARE 
IN  BAD  SHAPE.   RETURN  TO  FIRST  ORDER  METHOD  AND  TRY  AGAIN.   TF 
COURSE,  IF  NC  =  1  ALREADY,  THEN  THERE  IS  NC  HOPE  ANC  WE  EXIT  WITH 
KFLAG  =  -4. 


530  IF  (NQ.EC.l  )  GO  TO  540 


LCA 

3870 

LDA 

3E80 

LDA 

3890 

LCA 

3900 

LDA 

3910 

LDA 

3920 

LCA 

3930 

LDA 

3940 

LDA 

3950 

LCA 

3960 

LOA 

3970 

LCA 

3980 

LCA 

3990 

LDA 

4000 

LCA 

4010 

LCA 

4020 

LOA 

4030 

LDA 

4C40 

LDA 

4050 

LDA 

4060 

LCA 

4070 

LDA 

4080 

LCA 

4090 

LDA 

4  100 

LCA 

4110 

LCA 

4120 

LDA 

4130 

LDA 

4  140 

I  DA 

4150 

LDA 

4160 

LCA 

4170 

LOA 

4  180 

LDA 

4190 

LCA 

4200 

LDA 

4210 

LCA 

4220 

LDA 

4230 

LDA 

4240 

LCA 

4250 

LCA 

4260 

LCA 

4270 

LDA 

4280 

LDA 
■LDA 

4  290 
4  300 

LDA 

4310 

LDA 

4  3  20 

LCA 

4330 

LCA 

4340 

LDA 
■LUA 

4350 
4360 

LC«S 

4370 

LCA 

4380 

LCA 

4390 

LCA 

4400 

LDA 

4410 

I  DA 

4420 

I  DA 

4  4  30 

LCA 

4  4  40 

LCA 

4  4  50 

LCA 

4460 

LCA 

4470 

LCA 

4480 

L  DA 

4490 

LDA 

4  500 

LCA 

451C 

LDA 

4520 

LDA 

4530 

LDA 

4540 

■LDA 

*550 

LDA 

4560 

LDA 

4570 

LDA 

4  5  80 

LCA 

4  590 

•LDA 

4600 

LDA 

4610 
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NC  =  1  IDA  4620 

ID0U9  =  1  LCA  4630 

ASSIGN  570  TC  IRET  LDA  4640 

GC  TO  140  LCA  4650 

540  NCCLD  =  1  LOA  4660 

KFLAG  =  -4  LDA  467C 

GO  TO  220  LDA  468C 

550  KFLAG  =  -1  t CA  4690 

GC  TC  170  LDA  4700 

c      LDA  4  710 

C               THIS    SECTION    RESTORES    THE    SAVED    VALUES    OF    Y    AND    YL,     SCALING    ThE  LCA    4720 
C                Y    DERIVATIVES    AS    NECESSARY,     AND    THEN    PEiURNS    TO    THE    PREDICTS     LTOPI.DA    4730 

C LCA    4740 

560    H    =    HOLC*RACUM  LCA    4750 

H    =    AMAXUHP  IN,  AMINKH.HMAX)  I  LDA    4760 

570  PACUM  =  H/HCLC  LCA  4770 

Rl  =  1.  LDA  4780 

C  LCA  4790 

CO  580  J  =  2,K  LCA  4800 

Rl  =  R1*PACLM  LCA  4810 

C  LDA  4820 

DC  580  1  =  1.  NY  LDA  4830 

580  Y(J,II  =  SAVE(J,I)*R1  LDA  4840 

C  LDA  4e50 

C  LDA  4860 

CO  590  1  =  1,  NY  LDA  4870 

590  Yd, I)  =  SAVFd.I)  IDA  4880 

C  LDA  4890 

CALL  COPYZ  (YL.YLSV.LCCPYL)  LLA  4900 

IkEVAL  =  1  LDA  4910 

GO  TO  20C  LDA  4920 

6CC  KFLAG  =  -5  LDA  4930 

GC  TO  16C  LDA  494C 

c      LDA  4950 

C      THIS  SECTION  SCALES  THE  Y  DERIVATIVES  BY  R**J  LCA  4960 

c      LDA  4970 

610  Rl  =  1.  LCA  4980 

C  LDA  4990 

CC  620  J=2,K  LDA  5000 

PI  =  R1*R  LCA  5010 

C  LDA  5020 

CO  620  1=1, NY  LCA  5030 

620  Y(J,I)  =  Y(J,I»*R1  LOA  5040 

C  LDA  5050 

CC  TO  IRET,  (190,510)  LCA  5060 

c      LDA  5Q70 

C      THIS  SECTION  ALLJWS  FOR  RESTARTS  AFTER  rOLVINC  ANOTHER  PROBLEM  ORLDA  5080 

C      HAVING  TERMINATED  THE  CURRENT  COMPUTER  kUN.   SUBROUTINE  LDASAV  LDA  5090 

C      SAVES  THE  NECESSARY  VALUES  WHICH  ARE  INTERNAL  TO  LCASLB.   FOR  LDA  5100 

C      CCUBLE  PRECISION,  WITH  COPYZ  IN  SINGLE  PRECISION.  THE  NUMBER  OF  LCA  5110 

C      LCCATIONS  TC  PE  SAVED  &.ND  RESTORED,  LCOPYS  AND  LlOPYR.  MUST  fcE  LDA  5120 

C      SET  TC  56.  LDA  5130 

C      IT  IS  ASSUMED  THAT  IN  ADDITION  TO  THE  VARIAPLES  TN  THE  ARRAY  A  LCA  5140 

C      SAVED  BY  CALLING  LDASAV,  THE  USER  ALSO  SAVES  THE  ARRAYS  SAVE,  LCA  5150 

C      YLSV,  YMAX,  ESV,  AND  PW.  LDA  5160 

C  LCA  5170 
C      TC  RESTART  THE  USER  FIRST  CALLS  LDARST  TC  RESTORE  THE  VALUFi  SAVECLDA  5180 

C      eY  LDASAV,  THEN  RE-ENTERS  LDASUB  KITH  JSTAPT  <  C,  AND  WITH  THE  LDA  5190 

C      CTHER  PARAMETERS  THE  SAME  AS  RETURNED  FROM  THE  LAST  ENTRY  TJ  LDA  5200 

C      LCASUB,  PARTICULARLY  THOSE  ARRAYS  MENTIONED  ABOVE.  LCA  5210 

C      LCA  5  220 

ENTRY  LCASAWSAV)  LDA  5230 

LCOPYS  =  29  LDA  5240 

CALL  COPYZ  (SAV, A, LCOPYS)  LDA  5250 

CALL  COPYZ  (SAVE.Y.LCOPYY)  LDA  526C 

CALL  COPYZ  (YLSV, YL,LCOPYL)  LCA  5270 

RETURN  LCA  5280 

C  LCA  5290 

ENTPY  LCARST(SAV)  LCA  5300 

LCCPYP  =  29  LOA  5310 

CALL  COPYZ  (A,SAV,LCOPYR)  LDA  5320 

RETURN  LDA  5330 

C  LDA  5340 

C      LDA  5  350 

C  LCA  5360 
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c 

LDA 

5^70 

1  FC*MAT  (2I5,I2,1°2E10.2,7E14.6/(32X,7E14.6) > 

LCA 

5  2  80 

2  FCRMAT  (22X  ,1F7E14.6) 

LDA 

5  2  90 

3  FCPMAT  <«1   N  =',13,'    NL  =',13,'    KMSEPS  =  ',IPE9.2,'    T.NC  =  ' 

LCA 

5400 

i   ,£9.2,'    1-  =«,E9.2//) 
4  FCRMAT  («    NS    NW  Q     H',8X,'T    ',8X,'YU,*)  AMD  YL(*)'//) 

LDA 

5410 

LCA 

5420 

PNC 

LCA 

5430 

SLBfv3UTIN3  CCPYZ(S,Y,L) 
OIMENSICN  S(l)  ,Y(ll 

CCP 

10 

c 

COP 

20 
30 

c 

r.cp 

40 

c 

THIS  SUBROUTINE  CCFIE5  THE  ARRAY  Y,  OP  LENGTH  L,  INTO  TH!=  InFJV  S 

CC«= 

50 

I 

CCF 

60 
70 

IF(L.LE.C  )  RETURN 

C"P 

80 

CC  100  J=1,L 
IOC  S(J)  =  Y(J) 

CCF 

90 

COP 

100 

RETURN 

COP 

110 

ENC 

CCF 

120 

5LBR0UTINE  CERVAL  ( Y  ,  YL , T , N, NY, W , KERET ) 

DER 

10 

c 

DER 

20 

c 

THIS  SUBROUTINE  CALCULATES  THE  INTIAL  VALUES  uf    THE  DFRTVAT.VES 

DER 

30 

i 

IN  THE  GENERAL  CASE.   IT  IS  WRITTEN  SO  THAT  IT  SHOULD  WORK  IF  THF 

CEP 

40 

FIRST  NY  ECLA7I0NS  ALL  INVOLVE  DERIVATIVES.  l~    ATTEMPTS  TV    SOLVE 

CEP 

50 

c 

THE  FIRST  NY  EQUATIGNS  USING  NEWTON'S  METHCu,  BUT  SINCE  IT  7HES 

CER 

60 

c 

TC  EVALUATE  CF/DY'  BY  CALLING  JACMAT  IN  SUCH  A  WAY  a S  TO  MAKL  THE 

DER 

70 

c 

CF/DY  TERM  INSIGNIFICANT,  IT  IS  POSSIBLE  THAT  IT  MAY  FSR  FuP  THA 

TDEP 

80 

c 

REASON.   IT  MAY  FAIL  FOR  OTHER  REASONS,  AS  WELL.   IF  IT  00E5  FAIL 

CtP 

90 

c 

THE  USER  CAN  SUPPLY  HIS  uWN  VFRSION  OF  DLPVAL,  OP  MODIFY  TH4S 

CEP 

100 

c 

ROUTINE  IN  SUITABLE  FASHION.   THIS  ROUTINF  ASSUMES  THAT  VALuFS  CF 

CEP 

110 

c 

THE  LINEAR  VARIABLES  HAVE  BEEN  SUPPLIED  PREVIOUSLY.   IF  7HCSE 

CFR 

120 

c 

MUST  B^  SCLVEC  FOP  SIMULTANEOUSLY  WITH  THE  DERIVATIVES,  THP  USER 

DER 

13U 

c 

MLST  SUPPLY  HIS  0*N  VERSION  OF  DERVAL. 

DrR 

140 

c 

OER 

150 

c 

THE  CALLING  SEQUENCE  FIR  THIS  SUBROUTINE  IS 

DER 

160 

c 

DcP 

170 

c 

CALL  DcPVALlY,YL,T,N,NY~,W,KrPET) 

CER 

180 

s 

DtR 

190 

WHERE  THE  PARAMETERS  ARE  DEFINED  AS  FCLLGwS 

CER 

200 

c 

DER 

210 

c 

Y        -  SAME  AS  IN  LDASUB  ANO  SDESOL.   Yd, II  CONTAINS  THE 

CER 

220 

c 

INITIAL  VALUES  OF  THE  DEPFImDENT  VAPIA3LES.    THf 

DEP 

230 

c 

VALUES  OF  THE  DERIVATIVES  ARE  RETURNED  IN  Y(2,i). 

CEP 

240 

c 

YL        -  SAME  AS  IN  LDASUB  AND  SOEsGL.   THE  INITIAL  VALjES  Op 

CEP 

250 

c 

THE  LINEAR  VARIABLES  MUST  BE  SUPPLIED  TO  THIS  VFRSIONLEi 

260 

c 

T         -  INITIAL  TIME 

DEP 

27C 

c 

N        -  SAME  AS  IN  LDASUB,  TOTAL  NUMBER  OF  VARIABLES 

CER 

280 

c 

NY       -  SAME  AS  IN  LDASUB,  NUMBER  OF  CIFFFRENTIAL  EQUATIONS 

DEP 

290 

c 

AND  NONLINEAR  VARIABLES 

DEP 

20C 

c 

W        -  SCRATCH  ARRAY  W  FROM  THE  CALLING  SEQUENCE  OF  SDCSCL. 

CtP 

210 

c 

THIS  CAN  BE  USED  AS  NEEDED  IN  THIS  SUBROUTINE. 

DCP 

220 

c 

KERET    -  RETUPN  INDCATOP 

DER 

3  20 

c 

=0   NORMAL  RETURN 

DER 

240 

c 

=1   ERROR  RETURN 

DER 

2  50 

c 
c 

DER 

2  60 
370 
28C 

DIMENSION  Y(7,l),  YL(1),  W(l) 

DER 

c 

CEP 

390 

CG  100  1=1, NY 

DER 

400 

W(2*N+II  =  AMAXK  ABS(Y(  1,1  1  1,1.) 

Dc* 

410 

ICO  Y(3,I)  =  0. 

DEP 

420 

c 

DER 

430 

HINV  =  16.**20 

DER 

440 

KERET  =  C 

DEP 

450 

EFS2  »  NY/1.E8 

DER 

4  60 

EPS  =  SCPT(EPS2I 

DER 

4  70 

c 

DEP 

480 

DC  140  IT=1,10 

DER 

490 

c 

DEP 

500 

DC  110  1  =  1, NY 

OcR 

510 

110  Y(2,I  I  =  Y(2,I)/HINV 

DEP. 

520 
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C  DfcR  530 

CALL  DIFFUN  (  Y  ,  YL ,T , HI NV  ,  W  )  DEP  540 

CALL  JACMAT  IY,YL»T,HINV,-i.,NY,NY,5PS,W,W(N+l) ,W(3*N+1) )  DER  550 

NEWPW  -  1  DER  560 

C  DfcR  570 

CC  120  1=1, NY  DER  580 

120  W(  I  )  =  MI  )*HINV  DER  590 

C  DER  600 

CALL  NUITSL  <W( 3*N+1 ) ,W,W (N+l ) , NY, NY, E PS, W < 2+N+l) , NEWPW ,KRE 7 )  OER  610 

IF  (KPET.NF.O)  GO  TC  170  DER  620 

ER  =  0.  DfcF  630 

C  HER  640 

DC  130  1  =  1, NY  DFP  650 

Y(3,I)  =  Y(3,I)-W(N*I)  DER  660 

130  ER  =  ?R*U(N+I)/AMAX1(ABS(  Y(3,I  )  )  ,1.  )  )**2  HEP  670 

C  DER  680 

IF  (EP.LT.EPS2)  GO  TO  150  DER  690 

140  CCNTINUE  DEP  700 

C  DER  710 

GC  TC  170  DER  720 

C  DEP  730 

150  CC  160  1=1, NY  DEP  740 

160  Y(2,I  )  =  Y( 3,1)  DER  750 

C  DER  760 

RETURN  DER  770 

170  KERET  =  1  DER  780 

RETURN  DER  790 

ENC  DER  800 

SLEF3UTINE  JACMAT  ( Y, YL , T , HI NV, A2 ,N ,N Y , EPS , DY, F 1, PW)  JAC  10 

c      JAC  20 

C  JAC  3C 

C      SLBPOUTINE  JACMAT  IS  (USUALLY)  SUPPLIED  BY  THE  USER.   ITS  PbRPCSE  JAC  40 

C      IS  TO  EVALUATE  THE   J   MATRIX  NEEDED  WHEN  THE  CORRECTOR  EQUATION  JAC  50 

C      IS  SnLVEt  eY  NEKTON'S  METHOD.   THIS  VERSION  APPROXIMATES  JAC  60 

C      J   BY  NUMERICAL  CIFFERENCING  AND  USES  FULL  STORAGE  MOLE  JAC  70 

C       IN  AN  NXN  MATRIX.  JAC  30 

C  JAC  90 

C  JAC  100 

C      JACMAT  CALCLLATES  THE  MATRIX  JAC  110 

C  JAC  120 

C                 DF      A2  DF  JAC  130 

C         J   =    -----  —  jAC  140 

C                 CY      H   DY'  JAC  150 

C  JAC  160 

C      THE  CALLING  SEQUENCE  FOR  THIS  SUBROUTINE  IS  JAC  170 

C  JAC  180 

C      CALL  JACMAT(Y,YL,T,HINV,A2,EPS,N,NY,DY,F1,PW)  JAC  190 

C      WHERE  T|-E  PARAMETERS  ARE  DEFINED  AS  FOLLOWS.  JAC  200 

C  JAC  210 

C         Y        -  SAME  AS  IN  LDASUB  AND  IN  SDESOL.   ON  INPUT  TO  1'HIS  JAC  220 

C                    SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  Of  THE  JAC  230 

C                    DEPENDENT  VARIABLES  AND  THEIR  (SCALED)  DERIVATiVES.  JAC  240 

C         YL         SAME  AS  IN  LDASUB  AND  IN  SDESCl.   ON  INPUT  TO  7FIS  JAC  250 

C                    SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  OF  THE  JAC  260 

C                    LINEAR  VARIABLES.  JAC  270 

C         T        -  CURRENT  TIME  JAC  280 

C         HINV     -  1/H  ,  WHERE  H  IS  THE  CURRENT  STEPSIZd  JAC  290 

C         A2         A(2)  FROM  LDASUB.  JAC  300 

C         N        -  SAME  AS  IN  LDASUB,  TOTAL  NUMBER  OF  VARIABLES  JAC  310 

C         NY       -  SAME  AS  IN  LDASUB,  NUMBER  OF  DIFFERENTIAL   EQUATIONS  JAC  320 

C                    AND  NONLINEAR  VARIABLES  JAC  330 

C         EPS      -  L2  ERROR  CONSTANT  USED  IN  LDASUBl  JAC  340 

C         DY       -  ARRAY  OF  FUNCTION  VALUES  A'  CURRENT  VALUES  OF  THE  JAC  350 

C                    VARIABLES,  INPUT  TO  JACMAT.  JAC  360 

C         Fl       -  SCRATCH  ARRAY  OF  N  LOCATIONS  WHICH  CAN  BE  USED  EY  JAC  370 

C                    THIS  SUBROUTINE  IN  ANY  WAY  NEtCED.  JAC  330 

C          PW       -  J   MATRIX,  OR  APPROXIMATION,  CALCULATED  IN  JACMAT  ANDJAC  390 

C                    RETURNED  TO  CALLING  PROGRAM.   THIS  MATRIX  IS  USED  IN  JAC  400 

C                    SUBROUTINE  NUITSL  AND  STORAGE  MODE  MUST  AGREE  BETWEENJAC  410 

C                    THE  TWO  SUBROUTINES.  JAC  420 

C  JAC  430 

C      jAc  440 

DIMENSION  DYI1),  Y(7,l),  YL(1),  Fill),  °W(1)  JAC  450 
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NL  =  N-NY  JAC  460 

N'N  =  N*N  JAC  470 

C  JAC  460 

CC  100  1=1, NN  IAC  490 

100  Phd  )  =  0.  JAC  530 

C  JAC  510 

C  JAC  520 

CC  12C  J=1,NY  JAC  530 

F  =  Y(  1,J)  JAC  540 

E  =  Y(2,J)  JAC  550 

*  =  EPS*AMAXHEPS,ABS<F),ABS<E))  JAC  560 

Yd,  J)  =  Y(  1,J)+R  JAC  570 

Y(2tJ)  =  Y(2,J)-A2*R  JAC  580 

CALL  DIFFUN  (Y,YL  ,T,HINV,F1  )  JAC  590 

C  JAC  600 

CC  110  1=1, N  JAC  610 

110  PW(I  +  (J-lT*M  =  (  Fl(  I  )-0Y(  I  >  )/R  JAC  62C 

C  JAC  630 

Y(2,J)  =  F  JAC  64C 

12C  Yd, J)  =  F  JAC  650 

C  JAC  660 

IF  (NL.EC.O)  GO  Yd  150  JAC  670 

C  JAC  680 

CC  140  J=1,NL  JAC  690 

F  =  YL(JJ  JAC  700 

R  =  EPS**MAX1(EPS,ABS(F) )  JAC  710 

YL(J)  =  YUJJ+R  JAC  720 

CALL  DIFFUN  (  Y,  YL ,T , HI NV , F 1  )  JAC  730 

C  JAC  740 

CC  130  1=1, M  JAC  750 

130  PWU+U+NY-1  )*N)  =  (F1(II-JY(!)I/R  JAC  760 

C  JAC  770 

140  YL(J)  =  F  JAC  780 

C  JAC  790 

150  CCNTINUE  JAC  600 

RETURN  JAC  810 

PMC  JAC  820 

SUBROUTINE  NUITSL  ( PW, OY, Fl , N ,NY , EPS , YM AX, NEWPW , KRET )  NUI  10 

c      NUI  20 

C      THE  PURFCSE  OF  THIS  SUBROUTINE  IS  TO  SOLVE  A  NUI  30 

C      LINEAR  SYSTEM  OP  EQUATIONS  FOP  THE  NEWTON  ITERATES  WHEN  THF  NUI  40 

C      CORRECTGP  EQUATION  IS  BEING  SOLVED.   UPON  ENTRY  TO  THIS  SUBaCUTINENUI  50 

C      THE  SYSTEM  CF  EQUATIONS  TO  BE  SOLVED  IS    J  W  =  -c   ,  WFFRE  NUT  60 

C         J   IS  STORED  IN  PW  UPON  ENTRY  NUI  70 

C         WIS  RETURNED  IN  Fl  NU!  80 

C         -F  IS  STCRED  IN  DY  UPON  ENTRY  WUI  90 

C  NUI  100 

C      THIS  SUBROUTINE  IS  GENERALLY  SUPPLIED  BY  THF  USER,  ALTHOUGH  THERE  NUI  110 

C      ARE  SOME  STANCARQ  FORMS  AVAILABLE.   FOR  EXAMPLE,  THIS  \/ERSIJN  NUI  12C 

C       ASSUMES  THAT  PW  IS  STORED  IN  FULL  STOPAC-E  MODE  IN  AN  NXN  MATRIX.   NLI  130 

C      IF  NEWPW  =  1,  AN  LU  DECOMPOSITION  IS  DONE,  NEWPW  IS  SET  T0  ZERO    NUI  140 

C      Arc  FORWARC  AND  BACKWARD  SUBSTITUTION  FOR  THE  SOLUTION  IS  DGNE.    NUI  150 

C       IF  NEWPW  =  0,  ONLY  FORWARD  AND  BACKWARD  SUBSTITUTION  FOR  THE  NUI  160 

C      SOLUTION  IS  NECESSARY.  NUI  170 

C  NUI  180 

C      NOTE  THAT  THIS  VERSION  OF  NUITSL  REQUIPrS  THAT  PW  HAVE  N**2  +  2*fJ  NUI  190 

C       LOCATIONS  SINCE  2*N  LOCATIONS  ARfc  USCD  BY  THE  IMSL  LINEAR  EwL A TIONNUI  200 

C       SCLV5P.  NUI  210 

C  NUI  220 

C      NGTE  THAT  THE  PARAMETFRS  EPS  AND  YMAX  A^.E  LSEFUL  IF  aN  ITERATIVE   NUI  230 

C      METHOD  IS  USED  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS.  MUI  240 

C  NUI  250 

C      THE  CALLING  SEQUENCE  FOR  THIS  SUBROUTINE  IS  NUI  260 

C  NUI  270 

C      CALL  NUITSL (PW,DY  ,F1,N, NY, EPS, YMAX, NEWPW, KPET)  NUI  280 

C  NUI  290 

C       WHERE  THE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS.  NUI  300 

C  NUI  310 

C          CW          THE  J  MATRIX  CALCULATED  IN  SUBROUTINE  JACMAT  NUI  320 

C          DY        -  THE  RIGHT  HAND  SIDE  OF  THE  LINEAR  SYSTEM  TO  3E  SOLVEDNUI  330 

C          Fl        -  THE  SOLUTION  IS  RETURNED  IN  THE  ARRAY  Fl  NUI  340 

C          N         -  SAME  AS  IN  LDASUB,  TOTAL  NUMBER  OF  VARIAfiLtS  NUI  350 

C         NY         SAME  AS  IN  LDASUB,  NUM3FR  OF  C IFFERENT I A L  EQUATIONS   NUI  360 
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C  AND  NONLINEAR  VARIABLES 

C  EPS      -  L2  ERROR  CONSTANT  USED  IN  ..DASUB 

C  YMAX     -  MAXIMUM  VALUES  OF  Yd,  I  I  S  =  EN  UP  TO  THE  CURRENT  TIME 

C  NEWPW    -  INDICATES  WHETHER  A  NEW  J  MATRIX  HAS  BEEN  COMPUTED 

C  =1   INDICATES  A  NEW  J  MATRIX  HAS  BEEN  COMPUTED 

C  SINCE  THE  LAST  ENTRY  TO  NUITSL.   NEWPW 

C  SHOULD  BE  SET  TO  ZERO  IF  SOME  PREPROCESSING  NUI 

C  SUCH  AS  LU  DECOMPOSITION  MUST  BE  DONE  CN  A 

C  NEW   J   MATRIX. 

C  =0   INDICATES  THE  J  MATRIX  IS  THE  SAME  AS  V.HEN 

C  NUITSL  WAS  LAST  ENTERED 

C  KRET     -  RETURN  INDICATOR 

C  =0   NORMAL  RETURN 

C  =1   ERROR  RETURN.   SOLUTION  OF  EQUATIONS  CCULD 

C  NOT  BE  OBTAINED. 

C 

c      

DIMENSION  Pbd),  DY(1),  Fill),  YMAXd) 
NL  =  N-NY 

IF  (NEWPW. ECO)  GO  TO  100 
NEWPW  ■  0 
NN  =  N**2*l 
NNN  =  NN+N 

CALL  LUCATF  (PW.P W,N ,N ,0, Dl , D2, PWCNN ) , PW (NNN ) , Fl , I ER > 
IF  (  IER.EQ.O)  GO  TO  lOO 
KPET  *  1 
RETURN 
IOC  CALL  LUELMF  (PW,DY,PW( NN) ,N,N,F1  ) 
KRET  =  0 
RETURN 
END 

SUBROUTINE  CIFFUN <Y,YL,T ,HINV,DY ) 
C      

C  SUBROUTINE  CIFFUN  IS  SUPPLIED  BY  the  USER.   ITS  PURPOSE  IS  7C 

C  EVALUATE  THE  FUNCTIONS  AT  CURRENT  VALUES  OF  THE  VARIABLES. 

C 

C  THE  CALLING  SEQUENCE  FOR  THIS  SUBROUTINE  IS 

C 

C  CALL  DIFFUN(Y,YL,7,HINV,DY) 

C      WHERE  THE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS. 

C 

C  Y        -  SAME  AS  IN  LDASUB  AND  SDESOL.   ON  INPUT  TO  THIS 

C  SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  OF  THE 

C  DEPENDENT  VARIABLES  AND  THEIR  (SCALED)  DERIVATIVES. 

C  YL       -  SAME  AS  IN  LDASUB  AND  SDESOL.   ON  INPUT  TO  THIS 

C  SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  OF  THE 

C  LINEAR  VARIABLES. 

C  T        -  CURRENT  TIME 

C  HINV     -  1/H  .  WHERE  H  IS  THE  CURRENT  STEPSIZE 

C  DY       -  RETURNED  ARRAY  OF  FUNCTION  VALUES. 

C 

c      

DIMENSKN  Y(7,1),YL(1),DY(U 

C 

C      DEFINE  YOLR  FUNCTION  HERE 


C 


RETURN 
END 


NUI 

2  70 

NUI 

380 

NUI 

290 

NUI 

400 

NUI 

410 

NUI 

420 

NUI 

430 

NUI 

440 

NUI 

450 

NUI 

460 

NUI 

470 

NUI 

480 

NUI 

490 

NUI 

500 

NUI 

510 

NUI 

520 

•NUI 

530 

NUI 

540 

NUI 

550 

NUI 

560 

NUI 

570 

NUI 

580 

NUI 

5S0 

NUI 

600 

NUI 

610 

NUI 

620 

NUI 

630 

NUI 

640 
650 

NUI 

NUI 

660 

NUI 

670 
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Appendix  2 :   Examples 

Example  1:   This  example  is  the  problem  proposed  by  Gear  [3].   The 
system  of  equations  is 


7  4 

y±  -   S  +  (R-y,)Z  +  I     b4A   yA   =  0  ,  i  =  1,  2,  3,  4 


+  (R-y±)  +  I  b   y  =  0  ,  i  =  1 
j=l   J   J 

1  ? 

where  R  =  «•  2,  y  •   anc^ 

i=l  X 


l   4        ? 
S4J   (R-y,)   ,   and 
i=l 


y5  +  yl  y6  +  yl  y6  =  ° 

2y6  +  ye  -  H  +  vi  -  1  -  e_t  -  ° 
vi  -  v2  +  yi  n  ■  ° 

Vl  +  V2  +  5yl  y2  =  °  * 

In  the  above       b-n    =  b__   =  b„_   =  b,,    =  447.501 
11  22  33  44 

b12  =  "b34  "  b21  =  "b43  "  "  452-499 

b13  =  "b24  '  b31  "  "b42  =  "  47-499 

b14  -  "b23  '  b41  '  "b32  '  "  52-501   • 
The  initial   conditions   are 

y .  =  1  ,  i=l,  2,  3,  4  . 

y5  =  y0  =  i 

Vx  =  -  2  ,  V2  =  -  3  . 

3F 
Note  that  a  different  version  of  DERVAL  is  necessary  since  —  is  singular, 

3y 
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1C 


CIM 
CM 
DAT 

1  47 

2  44 
OAT 

1  1. 
CAL 
CAL 
CO 
I  G(  I 
DC 
Yd 
Yd 
Yd 
YL( 
YL( 
JSK 
CAL 
Pfil 
FCR 
STC 
END 


ENSI 
MSN 
A  GI 
.49? 
7.50 
A  \, 

e-4, 


(7,6),YL(2>  ,W<150>,3!  (16) 
/G(16) 

.501.- 452. 49 9, -4 7. 499, -52. 5 01, -4 52. 499, 44 7. 5  J 1,32.501, 
. 499, 5 2. 501, 4t7. 501, 452. 499, -52.  501,47.499,452. 499, 

L,N,REP3,HMAX,HMIN,H,T,TEND/8,fe,2,6,l.E-3,=.E2,l.E-l?, 


EC 


8  1  = 
)  = 
10  I 
,1  ) 
,5) 

if. 

2)  = 

F  = 
L  SD 
NT  6 
MAT( 
P 


2K  Y 

/CfcT 

/447 

,-47 

1/ 

NY,N 

u.  1 1 

RSET 

PSET 

1,16 

GI  (I  ) 

=  1,4 

=  -1. 

=  1. 

=  1. 

-2. 

-2. 
0 

ESCL(YfYL,T,TEND,NY,NL,M,JS<F,6,l,H,H*t*,l-mx,:>EP5,W) 
,JSKT 
'C 


.F3/ 
(207,256,-1,1 ) 
(208,256,-1,1) 


JSKF  =',I4) 


10 


20  S 


SIB 
CC* 
CAT 
DIN 
IF( 
TNT 
TCL 
CCN 
3  = 
s  = 

DC 
<;  = 


25 

30 


5C 


ICO 


ROUT 
M3N 
A  TO 
E.M5I 
T.EO 
ERM 
D  = 
TINU 

(Y( 

0. 


20 

c 


DY( 

CC 

DY( 

CCN 

CY( 

DY( 

DY( 

CY( 

RET 

fcNC 


30  I 
I  I  = 
25  J 
I  )  = 
TINU 

5)  = 

6)  = 

7)  - 

8)  = 
URN 


A'E  CIFFU.N,  (Y.YL.T.hlNV.DY) 

'CAT/G(4,4) 

.C/-13.459/ 

*  Y(7,l), YL(  l),DY( 1  I 

,TCLD)GO  Tl  10 

;   E  X  P  <  -T  ) 

id)  -»  Yd, 2)  +  Yd, 3)  +  Yd, 4)  )/2. 

1,4 

(R  -  Yd,  I)  )**2/2. 

1,4 

hINV*Y(2,I)  -  S  +  (R  -  Y(l,ill**2 

1,4 

CY(  I)  «■  G(I,J)*Y(1, J) 

•-INV*(Y(2,5)  +  Yd,  1)*Y(  2,6)  +  Y  (  2  ,  1 )  *Y<  1 ,  6  )  ) 
2.*Y(1,6)  +  Yd,6)**3  -  Y(l,l)  +  YL(l)  -  l.-TVTrR' 
YL(1)  -  YL(2)  +  Y(l,  1)*Y( 1,6) 
YL(1)  +  YL(2)  +  5.*YI1,1I  *  Yd, 2) 


SIB  ROUTINE  CERVAL(Y,YL,T,N,NY,rt  ,KECET) 

CIMEMSI  IN  Y(7,l), YL(  1) ,rt( 1) 

KERET  =  0 

CG  50  I=1,NY 

Y  (2  1 1  )  =  0. 

HINV  =  1. 

CALL  CIFFUN(Y,YL, T,HINV,WI 

CC     100  1=1, NY 

Y(2,I)  =  -MI  ) 

RETURN 

END 
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Example  2:   This  example  is  a  small  one,  contrived  to  illustrate  the 
possibility  of  derivatives  entering  in  a  nonlinear  fashion.   The  equations 
are 

y±  -   98Yl  +  98y2  =  0 

(yx)2  +  Y2  -  198y;L  +  e_t  y±  +   199y2  =  0 

v,  -  y,  -  y2  -  0 
The  initial  conditions  are 

yl  =  y2  =  1    »  Vl  =  2  ' 
Note  that  we  have  supplied  the  explicit  expression  for  the  Jacobian. 

Either  JACMAT  or  a  modified  version  of  DERVAL  must  be  supplied  as  the 

numerical  difference  approximation  to  the  Jacobian  causes  DERVAL  to 

fail. 
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DIMENSICN    Y(7,2) , YL ( 1) ,W( 50 » 

DATA    NtNY»NL»^tREPSiHMAXtHMIN»HtTiTEND/3i 2,  1, 2» 1* E-5, 25. » l.E-10, 
1    l.E-4,C..5C./ 
Y(l,l)     =    1. 
Yd, 2)     =    1. 
YL(1)     =    2. 
JSKF    =    0 

CALL  SDESCL  (YtYL,T tTEND ,NY,NL t M, JSKFf 6, 1 , hf HMlNt HMAX, REPS t W ) 
PRINT  6,JSKF 
6  FCPMAT( '0  JSKF  =  '  ,14) 
STCP 
ENC 

SUBROUTINE    CIFFUN    ( Y , YL , T , HINV, OY ) 
DIPENSICN    Y  (7,1), YL (  1) ,DY( 1) 
DATA    TCLC/-79.03/ 
IFU.EQ.TOLUGO    TO    10 
TKTERM    =    EXP(-T) 
TCLD    =    T 
10    CONTINUE 

DY(1)  =  Y(2,1)*HINV  -  98.*Y(1.1»  +  99.*Y(1,2) 

DY(2)  =  (Y(2,1)*HINV)**2  ♦  Y(2,2)*HINV  -  (196.  -  TMTE RM  )*Y ( 1 , 1 )  + 
1  199.*Y(1,2) 
DY (3)  =    YL(  1)  -  Y  (1,1)  -  Y( 1,2) 
RETURN 
ENC 

SUBROUTINE  JACMAT(YtYL,T,HINV,A2,N.NY,E°5,CY,Fl,PW  I 
CIMENSTCN  Y(7,ll, YL( 1),F1(1  ),DY(l),PW(Mt i) 
AH  =  A2*FINV 
CO  100  1=1, N 
CC  100  J=1,N 

ico  pmi ,j)   =  6. 

PV«  ( 1 1  1 )  =    -AH   -    98. 

Pin  ( 1  »2  )  =    9S. 

PW(2,1)  =    -2.*AH*Y(2,1)*HINV    -    198.    +    cXP(-T) 

Pk'(2,2)  =    -AH    +    199. 

IF(NL.LE.O)PETURN 

PM2.1)  =    -1. 

PM3,2)  =    -1. 

PM3.3)  =    1. 

RETURN 

ENC 
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Example  3:   This  is  another  contrived  example,  this  one  to  illustrate 
the  use  of  one  type  of  sparse  matrix  storage,  along  with  the  use  of 
iteration  to  solve  the  equations  (2.4).   The  system  of  equations  is 


A  y  +  B  y  =  0  , 


where 


A  = 


7 

0 

0 

-3 

0 

-1 

2 

8 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

3 

0 

0 

5 

0 

0 

0 

-1 

0 

0 

A 

0 

0 

0 

0 

-2 

0 

6 

B  - 


.3 

0 

0 

.1 

0 

-.2 

1 

3 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

10 

0 

0 

20 

0 

0 

0 

5 

0 

0 

6 

0 

lo 

0 

0 

57 

0 

100 

The  initial  conditions  are 
y1  =  1000 

y3  =  -25 

y4  -  io 
y5  =  o 

y6  =  -1000 
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The  matrix  storage  scheme  used  for  A,  B,   and  the  Jacobian, 
since  it  has  nonzero  elements  in  the  same  positions  as  A  and  B  ,  is 
that  outlined  by  Gustavson  [4],   Briefly,  one  stores  a  pointer  array 
(here  called  JS)  which  indicates  the  initial  position  of  new  elements 
in  two  other  arrays,  one  of  which  (here  called  JN)  gives  the  column 
number  of  the  element  stored  in  the  corresponding  position  of  the 
coefficient  arrays  (here  called  A  and  B) . 

Thus,  for  the  above  problem  the  arrays  stored  are 

JS:    1    4,   6-^7-^  9^~JL1       13 

4    6^2   """l  ±1 
A  :    7-3-1    8    2    1    5-3    4-1 


4 
-2 


B  : 


.3   .1  -.2 


20   10 


5  100   57 


The  elements  of  the  i —  row  are  stored  beginning  at  location 
JS(i)   of  the  array  A  ,  and  in  particular  A(JS(i)  +  k)   is  the  element 
in  the  i —  row  and  JN(JS(i)  +  k)  —  column  of  A  ,  for  k  =  0,  1,  ..., 
JS(i+l)  -  JS(i)  -  1  .   For  our  purposes  it  is  necessary  to  access  the 
diagonal  element  easily,  so  we  have  required  that  the  diagonal  element 
be  the  first  element  stored  for  a  given  row.   This  means  that 
JN(JS(i))  =  i  ,  i=l,  ...,  n  .   Note  that  JS(i)   is  the  number  of  nonzero 
elements  in  rows  1  through  i  -  1  , and  that  JS(n+l)  must  be  defined  as 
the  total  number  of  nonzero  elements. 

Problems  similar  to  the  above  arise  when  the  finite  element  method 
is  used  to  discretize  the  space  domain  for  time  dependent  partial  differ- 
ential equations.   Simple  modifications  to  the  subroutines  given  below 
should  permit  solution  of  large  problems  arising  in  that  fashion.   We 
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note,  however,  that  it  is  not  convenient  to  store  synmetric  matrices 
in  this  form  unless  all  nonzero  elements  are  stored.   Storage  of  only 
the  elements  of  the  lower  triangular  matrix  requires  one  to  reference 
columns  of  the  matrix,  which  are  not  readily  accessible.   Even  if  the 
entire  matrix  is  stored,  total  storage  requirements  for  matrices  arising 
in  finite  element  applications  is  still  considerably  less  with  this 
scheme  than  that  required  by  symmetric  band  storage  mode  [5], 
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ICO 

110 
120 


DI 

c: 

IN 

CA 

CA 

CA 

DA 

CA 

DA 

\ 

NF 

JE 

DC 

Y( 

DC 

A( 

e( 

JN 
DC 
JS 

PR 
PR 
pp 
PR 
CA 
PR 
ST 
CC 
FO 
FG 
FC 
FC 
FN 


NENSIC*  Y(7,6),M126),Yl(6),AD(12),rir,(12),JSC(7),JNC(12) 

MM3N  /CATA/A(12),B<12),N,JS17),JN(12> 

TEGER+2  JS.JN 

TA  T,TEi\D,F,JSKF  /  0  .  ,2  50  .  ,  1 .  E-5,  0  / 

TA  JSC/1,4,6,7,<;,11,13/ 

TA  JNC/1,3,5,2, 1,3,4, 1, 5,2,6, 4/ 

TA  AP/7.,-3.,-l.,8.,2.,l.,5.,-3.,4.,-l.,6.,-2./ 

TA  BC/.2,.l,-.2,3.,l.,l.,2T.,10.,o.,5.,lC0.,57./ 

TA  YT/1000.,0.,-2  5.,10. , J., -1000./ 


1  =  l\ 

=  JS 

100 

1.1  ) 

110 

I)  = 

I)  = 

(I  )  = 

120 
(I  )  = 
INT  4 
TNT  5 
INT  7 
INT  6 
LL  SO 
INT  3 
CP 

GKATt 
RMAT( 
R"AT( 
SMAT( 
RMAT( 


-  1 
) 


*  1 

C(NP1) 

I  =  1.N 

=  Y  I  (  I 

1=1, JE 

AC(I  ) 

ecn  ) 

JNCd 
I  =  1 ,  N  P 

JSC(T 
,N, (JS 
,(J|\(I 
,(A(I  ) 

,<e(  i) 

ESOL  (Y 
,JSKF 

•CRETLCN  FROM  SCESOL  vnI I  T H  JSKF  = • ,1 4 ) 

'  FJR  THIS  CASE  \=,II3//'  THE  JS  ARP.  «Y '  //(  1 2T 10 )) 

/•OTFE  J.N  APRAY«//(12I10)  ) 

/'OTFE  3  ARRAY'//(  12F10.2)  ) 

/■OTFE  A  ARRAY'//( 12F10.2) ) 


) 
1 
) 

(I  ) 
),I 
,1  = 
,1  = 
,YL 


,1=1, NP1 ) 
=  1 ,  J  E  ) 
If  JE) 

1  ,JE) 

,  T,TEN!D,.N,0»\t  JSKF,6,  l,H,l.E-h,U5.,l.E-4,<0 


300 
400 


SI 

CC 
IN 
CI 
CC 
CY 
JE 
JE 
DC 
CY 
CC 
CC 
R  E 
EN 


BrtZiUT 

TECER 

MENSI 

400 

(I )  = 

=    JS 

■-   JS 

300 

(I  )    = 

NTINU 
MTINL 

TURN 
C 


INC 
/CAT 
*2  J 
CN  Y 
1  =  1. 
C. 

( I ) 
(1  +  1 

J=JE 
CY( 

c 


CIFCJN     ( Y,YL,7,HINV,DY) 

A/M  12)  ,B(  12  ),N,  JS(7),  Ji\M  12) 

S,  JN 

(7,1), YL( 1) ,DY(  1) 

N 


)    -    1 

»JE 

I)     +    Y<2, JN( J) )*A(J)*HInV    +    ?(  J)*Y(  1,  JN<  J)  ) 


.     (PW,DY,F1,N,NY,EP'S,YVHX,N?WPW,KR 
CCMMDN    /CAT«/A( 12)  ,B(12),ND, JS(7) ,JN( 12 ) 


SUBROUTINE  NUITSL 


INTEGER +2    JS.JN 

DIMENSION  P«(l)  ,CY(  1)  ,F1(  1  ),  Yf'AX(  1) 

CATA  CMEC-»C*EGM  /1. 05,  .05/ 

KRET  =  0 

EFSS  =  EFS**2 

EFSA2    =    EP5S*.0001 

NCIT    =    N 

CC     100     1=1. NY 
100    Fid  )    =    CY  t  I)/Ow(  JS  1 1  )  ) 

CC  300  IT=1,NCIT 

RCh  =  0. 

CF  =  0. 

CC  200  1=1, NY 

JS  =  JS(I)  +  1 

JE  =  JS(I+1)  -  1 

FN:  =  DY(I) 

IF( JB.GT.JE  )GC  Tj  18C 

CC  150  J=JE,JE 
150  FN  =  FN  -  Pfc(  J)  *F1(  JN(  J)  ) 
130  FN  ---  FN/FW(je-l) 

FN  =  F\4C*EG  -  Fl  (I  )*GMEGM1 

flCH  =  Fill)  -  FN 
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O    =    CH    ♦    ( fiCH/Y^AXt I ) )**2 

RCH    =    »CH    +     ( ACH/AMAXiUBS (FN), EPS  I ) 
200    Fill)     =     FN 

IF(RCH.LT.EPSS)     RETURN 

IF(CH.lT.EPSA2)  RETURN 
300  CCNTINUE 

K9ET  =  1 

RETURN 

ENC 


*2 


ICO 


3  I  SPOUT  I NE  JACMAT(Y.YL,T,HIN\/,A2,N,NY,EPS,CY,F1,PW) 

CCMHCN    /CAT  A/41 12), 6(12 ) , NDUP, JS ( 7 ) , JN ( 12 > 

INTEGER*2  JS,JN 

CINENSICN  Y(7,l),YL(l),Fl(l),CY(l),PW(l) 

AK  *    -A2*HINV 

JE  =  JS1N+1  )  -  1 

CC  100  J=1,JE 

PW(J)  =  AH*MJ)  +  B(J) 

RETURN 

END 
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Example  4:   This  example  arises  from  a  nonlinear  reactor  dynamics 
problem  where  the  finite  element  method  is  used  to  discretize  the  space 
domain.   The  resulting  system  of  equations  has  the  form 

Ay-By  +  w(C  y)  y  =  0  , 


where  C  is  a  matrix  with  three  subscripts.   The  i —  equation  can  be 


expressed  as 


N       .    _  N    N  __ 

7   (a.   y.-b.   y.)+w  >    7  c.  ..  y .  y.  =  0  . 
J=l    J   J     J   J      J=l  k=l   J    J 


In  this  example  N  =  28  ,  and  there  are  at  most  seven  nonzero  elements 

per  row  in  A  and  B  .   The  nonlinear  term  y .  y,   appears  only  if 

J   K 

both  a. .   and  a,,   are  nonzero.   Therefore  a  different  type  of  sparse 
matrix  storage  is  used  for  this  problem. 

An  array,  K  ,  dimensioned  (28,  7)  is  used  to  store  (for  each  row), 
the  columns  subscripts  for  the  nonzero  elements.   For  convenience  in 
accessing  the  diagonal  element,  we  require  that  K(i,l)  =  i  .  We  can 
note  this  matrix  is  simply  the  connectivity  matrix  for  the  finite 
element  grid.   Then  the  nonzero  elements  of  A  and  B  ,  are  stored  in  the 
corresponding  portions  of  the  arrays  A  and  B  respectively.   If  there 
are  in  fact  less  than  seven  nonzero  coefficients  in  a  row,  the  remaining 
K(i,j)   are  set  to  zero. 

The  storage  for  C  is  somewhat  more  complicated.  C  is  symmetric 
(invariant  under  any  permutation  of  subscripts) .  The  nonlinear  term  of 
the  i —  equation  was  rewritten  as 
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N   N  _  N   N 

V  I       I      c        y,  yk  =  w  I       I      d    y  y   , 
j=l  k=l   1JtC  J   k     j=l  k=j   1Jk  J   k 


where 

k  <  j 

\jk   =   Pijk         k  =  J 

7ijk  +  Cikj  ■  k  >   J  ' 

The  coefficients  d     are  then  stored  in  a  (28,28)  array  C  in  the 
order  the  second  and  third  subscripts  are  given  here. 

(K(i,l),  K(i,l)) (K(i,l),  K(i,7)),  (K(i,2),  K(i,2)),  ...,  (K(i,2), 

K(i,7)),  ....  (K(i,7),  K(i,7))  . 

The  equations  can  then  be  written  in  the  form 


7  7   7 

^    [AiJ  ^(i,j)  "  Bij  yK(i,j)]  +  .1.    J.  Cim.k  ^(1,1)  yK(i,k)  "  ° 


where 


jk  2 


Because  of  the  large  amount  of  data  for  this  problem  the  input 
arrays  K,  A,  B,  and  C  are  simply  listed  along  with  the  programs  for 
this  example. 
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110 


120 


130 


135 
140 


c 
3 
10 
11 
12 
13 
14 
15 


CIMENSICN  Y<7.28)  ,WS<560> 

CC^MON  /CATA/A<28,7),B<28,7),C<2e,28),iN,NNZ,K<28,7> 

INTEGER*2  K 

CAT  A  TFNO tKHINtHM^X, EPS* ZOM EGA/. 1,1. E-12t.lt. 01 t 650. 903E-14/ 

CATA  NN',NY,NL/28,28,G/ 

NNZ  =  7 

CALL  ERRSET(2G7, 256, -1,11 

CALL  ERPSET(2C8,256,-i,l) 

N  =  NN 

¥    =  N 

NENC  =  NNZ*(NNZ  ♦  l)/2 

FRINT  10 

CC  110  1=1, NN 

READ  1,(K(I.J)  fJ=l,NNZ) 

FRINT  H,(K(I,J)  ,J=1,NNZ) 

"PINT  12 

CC  120  1=1. NN 

READ  2. (A(l ,J) ,J=1,NNZ) 

>( It  I  )  =  0. 

PRINT  15,(A(It J» » J=1,NNZ) 

Y(l,l)  =  1.E16 

PRINT  13 

CC  130  1=1, NN 

READ  2, (E( I, J),J=1,NNZ) 

PRINT  15 , ( e  < I , J) , J  =  1,NNZ) 

PRINT  \h 

CC  140  1  =  1,  NN 

READ  2,  (C(  I,J),J=1,NEND) 

CC  135  J=1.NEN0 

C(I,J)  =   C  ( I ,  J )  *ZGMEGA 

PRINT  15,(C(I,J), J=1,NEND) 

JSKF  =  0 

T  =  0. 

H  *  HI*IN*1000. 

CALL  SDESOL(Y,YL,T,TEND,NY,NL,K,  JCKF,  6,1,1-,  HMN.HMAX,   EPS. «<  ) 

PRINT  3, JSKF 

STCP 

FCRMAT(  1615) 

FGRMAT(7E11.4) 

FCPMATCO  JSKF  =  •,  13) 

FCRMATf  'IK  ARRAY' /) 

FCRMAT(8X,14I8) 

FCRMAT(///'OA(I, J )•/) 

FCRMAT<///«Oe<I ,J )'/) 

FCPMAT(///'OC(I,J ) '/ ) 

FCRMAT(8X,1F7E16.6) 

END 


SLBROUTINE  CIFFUN(Y,YL,T,HINV,DY) 

CCMMGN  /CAT A/ A (28, 7) ,6(28, 7),C( 28,28) ,N»NNZ,«(28, 7) 

INTEGER*:  K 

DIMFNSKN  Y(7,l)  ,  YL(1)  ,DY(  1) 

CO  400  1  =  1, N 

CY(I )  =  0. 

CO  300  J=1,NNZ 

IF(K(I,  JJ.LE.OJGC-  TO  310 

OY(I)  =  CY(I)  +  Y(2,K( I, J)  )*A(I , Jl+HINV 
1  C(I,J)+Y(1,K(I,J))*Y(1,K(!,1)) 
30C  CCNTINUE 
310  L  =  NNZ 

CC    360    J1=2,NNZ 

IF  (K(  I,  JD.LE.OJGC    TO    400 

DC    350    J.~  =  J1,NNZ 

L    =    L    +    1 

IF(K(  I, J2).LE.0)G0    TO    350 


BU,J)*Y(1,K<  I,  J)  )    + 


DY(I  )  =  CY( I ) 
350  CCNTINUE 
360  CCNTINUE 
400    CCNTINUE 

RETURN 

END 


+  C (I  ,L)*Y(l,Kt I, Jl)  )*Y(1,K(T  ,  J2) ) 


SLePOUTINE  JACMAT(YtYL,TtHINV,A2tN,NY,EPS,CY,Fl,PW) 
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96 
99 


ICO 


200 
3G0 


280 
281 


285 

287 
288 


COMMON  /CATA/A<28,7),B(28,7),C<28,23),NOfNNZ,K(28,7) 

INTEGERS  K,P 

DIMENSION  V(7tl).VLIll,Fl(l  )  ,DY <  1 )  , PW ( NY  ,  1  ) 

CIPENSI1N  P(7,7) 

CATA  P/49*0/,lNITP/0/ 

IF(INITP.EC.NNZ)GO  TQ  99 

IMTP  =  NNZ 

CO  98  1=1, NNZ 

CC  98  M  =  L,NNZ 

P(L,M)  *  ft    +  (L 

CCNTINUE 

AH  =  -A2*HNV 

CO  300  1  =  1, NY 

DG  300  J=1,NNZ 

PWU  ,J)  =  Ah*A(  I,  J) 

CC  100  L=1,J 
IF(K(I,U.LE.0)GO 

PU(I,J)  =  PMI.J)  +  C(T  ,P(L,  J))*Y(1,K(!,L)  ) 
CO  200  N=J,NNZ 
IF(K(I,y).LE.OJGC 
PfcU,  J) 


1)*(2*NNZ  -  L )/2 


-  b( :, j) 

TO  295 

+  C(T  ,P(L, J) )*Y(1,K(I 


PkU  ,ji 

CCNTINUE 
CCNTINUE 
RETURN 
END 


TO  295 

+  C(I  ,P(  J,Vt)  )*Y(1,K(I  ,M>) 


SUBROUTINE  NUITSl (Pd ,0Y ,c 1 , N, NY , EPS , YMAX , NEWPW, KRET ) 

DIMENSION  PMNY.l ) ,DY(1) , F I ( 1  ) , YMAX ( 1 ) 

CCPMCN  /CATA/A<28,7),B(28,7),C(23,28>,'J0,NNZ,M23,7) 

CATA   SPC,SPCI"l/l.05,.05/ 

INTEGEP*2  K 

KR£T  =  0 


INTEGEP*2 
KRE"""  =  0 
EFSS  =  EFS**2 
EPSA2  =  FPSS+.0001 
NCIT  =  N 
CC  281  1  =  1, NY 
Fl(I)  =  CY(I)/PW< 1,1) 
CC  287  17  =  1, NOTT 
RC(-  =  0. 
CI-  =  0. 
CO  285  1=1, NY 
FN  =  DY(I) 
CG  284  J=2,NNZ 

IF(K(  I,  J  J.LE.O.OP.KC,  J).GT. 
FN  =  CN  -  PMI,  J)*FHK<  !,  J)  ) 
284  CONTINUE 

FN  =  FN/FM!  ,1) 

FN  =  FN*SPb 

ACH  =  Fid  ) 

CH  =  CK 

PCH  =  R_ 

Fl  (I  )  =  FN 

IF(RCH.LT.EPSS)Gn 

IF(Ct-.LE.EPSA2)GQ 

CCNTINUE 

KRET  =  3 

CCNTINUE 

RETURN 

END 


,J).GT.NY)GO  TJ  284 


SPOyi*Fl< I) 
:1 (I )  -  FN 

*    *    ( ACH/YMAX( I ) )**2 
5CF  +  (ACH/AMAX1(A6S(FN 

»  FN 


),EPS) 1^*2 


TG 

TO 


288 
288 
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INPUT  CATA  FOR  EXAMPLE  4 


1 

6 

2 

2 

1 

6 

7 

8 

3 

3 

2 

8 

4 

4 

3 

£ 

9 

10 

5 

5 

4 

1C 

21 

6 

11 

7 

2 

1 

7 

2 

6 

11 

12 

13 

a 

8 

3 

2 

7 

13 

9 

4 

9 

4 

8 

13 

14 

15 

10 

10 

5 

4 

9 

15 

22 

21 

11 

6 

16 

12 

7 

12 

7 

11 

16 

17 

18 

13 

13 

8 

7 

12 

18 

14 

9 

14 

9 

13 

16 

19 

20 

15 

15 

10 

9 

14 

20 

23 

22 

16 

11 

28 

17 

12 

17 

12 

16 

28 

27 

18 

18 

13 

12 

17 

27 

26 

19 

19 

14 

18 

26 

25 

20 

20 

15 

14 

19 

25 

24 

23 

21 

5 

1C 

22 

22 

21 

10 

15 

23 

23 

22 

15 

20 

24 

24 

23 

20 

25 

25 

20 

19 

26 

24 

26 

19 

16 

27 

27 

18 

17 

28 

26 

28 

16 

27 

17 

14 


8.3776E  02  8.3776E  02  4.1888E  02  0.0         0.0         0.0         0.0 

5.0265E  03    4.1888E    02    2.0944E    03    2.5133E    03    2.0944F    03    4.1888E    02    0.0 

1.6755E  03    4.1888E    02    1.6755E    03    4.1888E    02    0.0 

5.0265E  03    4.1888E    02    2.0944E    03    2.5133E    03    2.0< 

1.4661E  03    4.1888E    02    1.4661E    03    3.14166    02    0.0  0.0  0.0 


5.0265E    03    4.1888E    02    2.0944E    03    2.5133E    03    2.0944E    03    4.1888E    02    0.0 

_    03    4.1888E    02    1.4661E    03    3.1' _ 
1.C891E    04    2.<322E    03    4.1888E    03    2.0944E    03    8.37765    02    0.0  0.0 


2.8484E    04  2.5133E    03  4.1888E  03    6.2832E    03    6.7021E  03  6.2832E    Oi    4.1883E    03 

2.1782E    04  1.6755E    03  2.0944E  03    4.1888E    03    5.8643E  03  4.16886    Oi    2.0944E    03 

2.6464E    04  2.5133E    03  4.1888E  03    6.2832E    03   6.7021E  03  ^>.2832E    03    4.18385   03 

1.9C59E    04  1.4661E    03  2.0944E  03    4.1888E    03    5.13135  03  3.1416E    03    1.5708E    03 

2.3457E    04  2.«322E    03  5.0265E  03    8.3776E    03    6.2832E  03  0.0                     0.0 

5.3616E    04  6.7C21E    03  8.3776E  03    1.0472E    04    1.C891E  04  1.0472E    0-*    8.3776E    03 

4.6914E    04  5.8642E    03  6.2832E  03    8.3776E    03    1.0053E  04  8.3776E    03    6.2832E    03 

5.3616E    04  6.7021E    03  8.3776E  03    1.0472E    04    1.C891E  04  1.0472E    04    6.3776E    03 

4.4663E    04  5.1313E    03  6.2832E  03    8.3776F    03    8.9535E  03  8.4823E    0-.    6.2046E    03 

3.2515E    04  5.0265E   03  5.1836E  03    1.0812E    04    1.0472E  04  0.0                     0.0 

5.62C8E    04  1.C891E    04  1.0812E  04    1.1958E    04    1.1958E  0*  1.0812E    04    0.0 

8.0582E    04  1.0052E    04  1.0472E  04    1.0812E    04    1.33125  04  1.3212F    04    2.12845    04 

5.62C8E    04  1.0891E    04  1.G312E  04    1.1958E    04    1.1958E  04  1.0612E    0*    0.0 

6.C750E    04  8.S535E    03  1.0472E  04    1.0812E    04    9.2481E  03  1.0348E    04    9.8764E    03 

3.7699E    03  3.1416E    02  1.5708E  03    1.8850E    03 

2.9531E    04  1.6650E    03  3.141&E  03    6.2046E    03    8.7179E  03 

6.1889E    04  e.7179E    03  8.4823E  03    9. 6764E    03    1.2468E  04 

6.33035    04  1.2466E    04  1.03V8E  04    8.8357E    03 

5.7962E    04  9.2481E    03  1.1958E  04    1.4726E    04    8986575  03 

5.4719E    04  1.1958E    04  1.3312E  04    1.7671E    04 

9.4719E    04  1.3312E    04  1.1958E  04    1.4726E    04    1.7671E  04 

5.3014E    04  5.1826E    03  1.4726E  04    1.1958E    04 

■1.5816E    09  1.1720E    09  1.0449E  09    0.0                     0.0  0.0                     0.0 

1.0449E    09    0.0 
0.0  0.0 

■3.9622E    09  1.0449E    09  6.3535E  08    4.4338E    09    6.3535E  08  1.0449E    09    0.0 

■4.3361E    09  1.0449E    09  1.8355E  09    1.4879E    09    0.0  0.0                      0.0 

■6.7925E    09  4.560SE    09  6.7778E  09    6.35355    08    1.1720E  09  0.0                      0.0 

1.5222E    10  4.4338E    09  6.7778E  09    1.90616    09    1.1212E  10  1.9061E    09    6.7776c    09 


-1.5816E    09    1.1720E    09    1.0449E    09    0.0  0.0 

-3.9622E    09    1.04495    09    6.3535E    08    4.4338E    0*    3.3535E    08 
-3.1631E    09    1.0449E    09    2.3440E    09    1.0449E    09  0.0 


-1.3585E  10  2.2440E  09  6.3535E  08    6.7778E    09  9.12185  09  6.77785    05    6.3535E  08 

-1.5223E  10  4.4338E  09  6.7778E  09    1.9061E    09  1.1212E  10  1.90615    09    6.7778E  09 

-2.4104E  10  1.8355E  09  6.3535E  08    6.7778E    09  7.3355E  09  8.44466    09-6.0318E  08 

-1.3995E  10  4.E609E  09  7.9498E  09    1.3556E    10  1.9061E  09  0.0                      0.0 

-2.9627E  10  1.1212E  10  1.3556E  10    3.1768E    09  1.7989E  10  3.17666    09    1.3556E  10 

-2.7989E  10  9.1218E  09  1.9061E  09    1.3556E     10  1.5900E  10  1.3556E    lu    1.9061E  09 

-2.9627E  10  1.1212E  10  1.3556E  10    3.1768E    09  1.7989E  10  3.1768E    09    1.3556E  10 

-4.8828E  10  7.3355E  09  1.90616  09    1.3556E     10  1.0212E  10  1.16215    10    2.04066  09 

-3.5210E  10  7.9498E  09  1.3692E  10    1.6043E     10  3.1766E  09  0.0                      0.0 

-8.2408E  10  1.7989E  10  1.6043E  10-3.6945E    08  2.40606  10  1.31036    1C    0.0 
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-7.2221E 

-8.24C3E 

-8.7119E 

-8.2637E 

-4.7C95E 

-1.01266 

-1.32805 

-1.2246E 

-1. 5263E 

-1.4877E 

-7.4644E 

4.1888E 

4.1388= 

0.0 

0.0 

2.5123E 
1.3S63E 
0.0 

2.2340E 
8.3776E 
1.39636 
2.7925E 
0.0 

2.5123E 
1.3963E 
0.0 

2.2240E 
4.18885 
1.3963E 
0.0 
0.0 

6.70215 
1.2566? 
5.5850= 
1.1170E 
1.6755E 
8.3776E 
1.11705 
2.2340= 
1.3404E 
5.56506 
0.0 

1.9548E 
1.6755E 
8.3776= 
1.1170E 
2.2340E 
6.7C21E 
2.7925E 
0.0 

1.9548E 
1.4242E 
1.1170E 
1.9548E 
2.0718E 
3.1835E 
5.5132E 
1.9548E 
2.9G95E 
2.8484E 
2.23406 
0.0 

4.46805 
3.1835E 
2.5133E 
1.9549t 
3.9095E 
1.42426 
1.1170E 
0.0 

4.46806 
1.3822E 
1.95486 
0.0 
4. 747  3E 


10  1 
10  1 
10  1 
09  I 
10 
11 


11 
11 
11 
11 


10    1 


02 
02 


03 
02 

C? 
02 
02 
02 

03 
02 

03 
02 
02 


5 
2 
0 
0 
0 
2 
0 
9 
2 
■> 

6 
o 
o 

2 
6 

9 
2 

0 
0 

02  2 

03  1 
02    0 


o: 

04 
02 
02 
C3 
04 
02 

03 

04 
02 
03 

03 
03 
02 

03 
04 
03 
03 
03 
04 


03    1 
03    0 

03  4 

04  0 

03  1 
0 

02  4 

04  5 

03  1 
03    0 

03  4 

04  0 
03    1 


03 
04 
03 

03 


.59CCE 

.7989E 

.021 2E 

.4879E 

.89535 

.40C5E 

.6822E 

.4796E 

.406CE 

.14765 

.2692? 

.5850E 

.7925E 

.C 

.0 

.0 

.7925E 

.0 

.7738E 

.7925E 

.79255 

.0 

.0 

.c 

.79255 

.0 

.77386 

.7925E 

.7925E 

.0 

.C 

.2340E 

.U70E 

.0 

.7925E 

.23406 

.5S50E 

.0 

.7925E 

.0 

.7925E 

.C 

.513  26 

.2340E 

.5650E 

.0 

.7925E 

.0 

.7925E 

.C 

.2566E 

.0 

.c 

.0 

.2  736E 

.565CE 

.2962E 

.0 

.46805 

.0 

.1170E 

.0 

.1388E 

.5850E 

.29  6  2E 

.C 

.46805 

.0 

.1170E 

.C 

.0944c 

.0 

.0 

.0 

.0492E 


10 

10 

10 

09- 

09 

08 

09 

10- 

10 

10 

10 

02 

02 


02 

02 
02 
02 


02 

32 

02 
02 


03 
J3 

02 
02 
02 

03 

02 

03 
03 
02 


02 
03 


03 
03 
03 

03 

03 

03 
03 
03 

03 

03 

J  3 


03 


3.1768E 

09 

1.3103E 

10 

3.1768E 

09 

6.0318E 

08 

8.4446E 

09 

1.1621F 

10 

2.4658E 

09 

3.6945E 

08 

1. 1476E 

10 

2.406JE 

10 

3.3929E 

09 

2. 7925E 

02 

0.0 

0.0 

0.0 

8.3776E 

02 

0.0 

0.0 

0.0 

1.U70E 

03 

0.0 

0.0 

0.0 

8.3776E 

0? 

0.0 

0.0 

0.0 

5.5850E 

02 

0.0 

CO 

0.0 

1.9548E 

03 

CO 

0.0 

0.0 

1.3963E 

03 

0.0 

CO 

1.39636 

03 

1.1170E 

03 

0.0 

0.0 

1.1170E 

03 

1.3963E 

03 

0.0 

0.0 

1.3963E 

Oj 

1.1170E 

03 

CO 

CO 

0.0 

2.9095E 

03 

CO 

0.0 

0.0 

2.0716E 

03 

0.0 

0.0 

2.2240E 

03 

2. 79 2  56 

03 

0.0 

CO 

1.9548E 

03 

3.0718E 

03 

CO 

CO 

2.2  340E 

03 

2.7925E 

03 

0.0 

0.0 

0.0 

CO 

CO 

0.0 

CO 

1.31035 

2.*060c 

1.60436 

2.89  5  3F 

2.04095 

1.3398E 

2.3750= 

7.2533E 
•2.5688E 

3.3929E 
-3.6945F 

0.0 

0.0 

0.0 

0.0 

J.O 

0.0 

5.5850E 

8.3776E 

2.7925F 

0.0 

0.0 

0.0 

0.0 

0.0 

5.58505 

3.37765 

0.0 

0.0 

0.0 

0.0 

3.37765 

0.0 

1.3963E 

1.U70F 

2.513  35 

0.0 

1.95485 

6.42285 

1.9549E 

0.0 

8.3776E 

0.0 

2.5132E 

0.0 

1.95436 

6.4228E 

1.95486 

0.0 

3.3776F 

0.0 

3.&303F 

1.1170E 

0.0 

2.7925F 

4.1868E 

0.3 

3.63035 

1.0612c 

3.6303E 

0.0 

2.5133E 

0.0 

4.18885 

0.0 

3.6303E 

1.0612E 

3.6303E 

0,0 

2.51335 

0.0 

0.0 

1.9548F 

0.0 

4.4680E 


10  1, 
10-3, 
1.0 

39 
0  9 
10 

in 

09    2.3750E    10 

09 

09- 

08 


.147bE  10 

!. 69455  08 

1.47986  10 

..4C055  08 

h. 68226  09 


1.14765    10    1.62606     10 
1.6C43L    1C    CO 
3.46585    09    1.33986    1C 


3.56885    09 


0 
0 
0 
0 
5 
0 
02  1 
02    1 


02 


0 

0 

1 

0 
5 
0 
02  1 
02  1 
0 
J 
0 
0 
0 
0 
6 
0 
0 
0 
2 
2 
4 
0 
1 
1 
0 
0 
2 
2 
2 
0 

02  1 
0 
-> 

0 

2 
0 
0 
0 


0  2 

03 
03 
03 

03 
02 

03 


Oj 

03 
02 
03 


03 
0  3 

02 
03 

02 
04 
03 

03 

03 

03 
04 
03 

03 


03 
Oi 


.0 

.0 

.0 

.0 

.58505    02 

.0 

.1170=    Oj 

.39636    02 

.C 

.0 

.3963F    02 

.0 

.585CE    02 

.0 

.1170F    03 

.3963=    02 

.0 

.0 

.0 

.0 

.0 

.0 

.98135    02 

.0 

.0 

.0 

.65295    03 

.6529E     03 

.46806    03 

.0 

.67556    03 

.67555    03 

.0 

.0 

.65296    03 

.65296    03 

.23406    03 

.0 

.6755=    03 

.0 

.51226    03 

.0 

.3510F    03 

.C 


0.0 
CO 
0.0 

0.0 

0.0 

2.7925- 

5.56505 

5.585C6 

0.0 

0.0 

0.0 

0.0 

CO 

2.7925c 

5.56506 

5.56505 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

2.22405 

0.0 

0.0 

2.5133  = 

1.9548E 


0 

1 

0 

0 

0 

02    9 

02    0 

02    0 

0 


02 


02 

C2 
02 


03 


02 

0  3 


.22845    03 

.32346    0? 

.81912    03 

.0 

.35iOE 

.25105 

.0 

.0 

.3284E    03 

.3234E    03 

.9C95E    03 

.0 

.3510=    03 

.0 

.ie886    03 

.0 

.37366    03 

.0 


1.9 
1.1 
l.« 
2.2 
0.0 
0.0 
2.5 
0.0 
9.2 
1.1 
0.0 
CO 
2.5 
1.3 
0.0 
2.9 
0.0 
0.0 
4.1 
2.6 
6.1 
1.9 
4.4 
3.9 
O.C 
CO 
4.1 
0.0 
2.5 
1.9 
0.0 
0.0 
4.1 
2.2 
0.0 


548^ 
1705 

5*86 
240E 


0; 
02 
03 
Oj 


6 

0 
0  2  1 
Oj    1 

1 
0 


1336    03 


7766 
17CE 


1325 
9625 


o: 

Oj 


0: 

Co 


J91 


88  8  6 
2026 
4265 
548E 
6d05 
0956 


5 
0 
6 
0 

c 

0 
2 
0 
0 

03  2 
3 
0 
Oj  1 
05  2 
Oj  2 
03  0 

02  5 

03  2 


8886  Cj 


1322 

5485 


88  86 
240c 


Oj 

0. 


Oj 

Oj 


.0 

.396: 

.0 

.0 

.0 

.77366  02 

.0 

.0 

.0 

.27766  02 

.0 

.0 

.0 

.7738E  02 

.0 

.0 

.0 

.1336=  02 

.0 

.0 

.0 

.6755E  03 

.0 

.0 

.29625  03 

.67555  02 

.0 

.3058=  03 

.3  7  76  5  0  2 

,9813t  02 

.0 

.81515  02 

.3962c  03 

.67556  03 

.0 

.3053E  03 

198136  02 

.0 

.0 

.0 

.09445  03 

•  u 

.0 

.0719E  0? 

.35106  01 

.0 

.03325  04 

.5132=  02 

.3726=  03 

.0 

.16625  03 

.U7185  03 

.351C6  03 

.0 

.03325  04 

.0 

.3736E  03 

.0 

.0 

.0 

.0 

.0 

.0 
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1.4242E 

04 

8.9361E 

8? 

4.7473E 

03    0.0 

0.0 

4.7473E 

03 

0.0 

4.1868E 

03 

2.234CE 

0.0 

0.0 

0.0 

0.0 

2.3736E 

03 

0.0 

0.0 

Q.O 
O.O 

0.0 

0.0 

0.0 

CO 

0.0 

0.0 

0.0 

2.3736E 

03    2.22 

40E 

03 

0.0 

2.7646E 

04 

0.0 

4.4680E 

03    0.0 

0.0 

0.0 

4.1888E 

03 

3.9095E 

03 

1.95485 

03 

0.0 

0.0 

0.0 

1.02 

32E 

04 

4.0492E 

03 

0.0 

0.0 

0.0 

0.0 

4.188eE 

03    2.3736E 

03    0.0 

0.0 

6.9813E 

03 

0.0 

0.0 

8:8 

0.0 

1.7872E 
0.0 

04 

1.4242E 

04 

e.9361E 

03 

4.7473E 

03    0.0 

4.7472E 

Oj 

4.1888E 

02 

2.234CE 

03 

0.0 

0.0 

0.0 

0.0 

2.3736E 

03 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

2.3736E 

03    2.22 

40E 

03 

CO 

1.3823E 

04 

O.C 

4.4680E 

03    0.0 

0.0 

0.0 

0.0 

1.9548E 

03 

1.9548E 

33 

0.0 

0.0 

0.0 

4.1888E 

03 

4.0492E 

03 

0.0 

g.o 

0.0 

0.0 

4.1888E 

03    2.3736E 

03    0.0 

0.0 

6.9813E 

03 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

O.Q 
0.0 

0.0 

0.0 

g.o 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 
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